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The  g»thematical  modeling  techniques 


to  predict  atmospheric 


dispersion  of  heavy  gases  are  reviewed  and  critiqued.  These  dispersion  processes 
are  characterized  by  three  ‘regimes  -etfr  phases^  fi)  . 

^^-^-nearyfield,  negative  buoyancy-dominated  flow  regime^  ' 


intermediatewfield,  stably  stratified  shear  flow  regime ^ 

>far±field,  passive  turbulent  diffusion  regime. 

Mathematical  models  of  laboratory  and  natural  gravity  currents  are  used  to  describe 
the  negative  buoyancy-dominated  regime  flow  and  dilution  processes.  A  laboratory 
investigation  of  instantaneous  heavy  gas  releases  in  calm  air  demonstrates  scaling  ^ 
from  small  release  volumes  (0.5  ffij^pEo  large  field  (2000  (PT)  field  releases,  and  is 
used  to  develop  a  mathematical  model  for  the  buoyancy-domiffated  regime.  Laboratory- 
data  from  stratified  shear  flow  mixing  experiments  are  used  to  develop  a  model  for 
vertical  diffusion  of  heavy  gases  in  the  atmospheric  surface  layer.  The  model 
concepts  are  consistent  with  the  limiting  passive  behavior  of  demonstrated  air 
pollution  models.  An  interactive  computer  model  for  heavy  gas  dispersion  (DEGADIS) 
is  developed.  The  DEGADIS  model  accounts  for  the  three  regimes  of  heavy  gas  dis¬ 
persion  processes  and  can  be  used  to  simulate  instantaneous,  steady  state,  and 
transient  releases.  The  model  provides  predictions  of  downwind  concentration  decay 
which  are  in  good  agreement  with  the  full  range  of  field  experimental,  data  currently 
available,,  and  is  recormtended  for  incorporation  in  the  Coast  Guard  Hazard ^Assessment 
J26mputer  ^aystem.^HACS) . 
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vertical  turbulent  diffusivity  (m  /s) 
kinetic  energy  (Nm) 
von  Karman's  constant,  0.35 
constant  in  Equation  11-47 
constant  in  Equation  11-48 

source  length  (m),  or  turbulence  length  scale  (m) 

1/3 

characteristic  length  V/ 

total  cloud  mass  (kg) 

total  mass  of  air  in  the  cloud  (kg) 

total  mass  of  contaminant  in  the  cloud  (kg) 

initial  cloud  mass  (kg) 

molecular  weight 

mass  rate  of  air  entrainment  into  the  cloud  (kg/s) 
mass  rate  of  water  transfer  to  the  cloud  from  the  water 
surface  under  the  source  (kg/s) 


number  of  observers 


Nu 


«s 

Q. 

^*max 

«s 

qo 

ql 

R 


^max 


Nusselt  number 
potential  energy  (Nm) 

Prandtl  number 
atmospheric  pressure  (atm) 
empirical  constant  in  Table  A.l 
vapor  pressure  of  water  (atm) 
cloud  momentum  (kg  m/s) 

momentum  of  head  in  density-driven  flow  (kg  m/s) 

momentum  of  tail  in  density-driven  flow  (kg  m/s) 

virtual  momentum  due  to  acceleration  reaction  (kg  m/s) 

volumetric  release  rate  (m  /s) 

2 

source  mass  flux  (kg/m  s) 
volumetric  entrainment  flux  (m/s) 

flux  of  ambient  fluid  into  front  of  gravity  current  head 
(m/s) 

rate  of  heat  transfer  from  the  surface  (J/s) 

p 

atmospheric  takeup  flux  (kg/m  s) 

maximum  atmospheric  takeup  flux  of  contaminant  (kg/m2  s) 

surface  heat  flux  (J/m  s) 

vertical  heat  flux  (0/m  s) 

empirical  constant  in  Table  A.l 

gas  source  cloud  radius  (m) 

inner  radius  of  head  in  density-driven  flow  (m) 

value  of  R  when  (it  R  Q*)  is  a  maximum  (m) 

maximum  radius  of  the  cloud  (m) 


Rp  primary  source  radius  (m) 

R*  nondimensionalized  radius  (R/Jl) 

Ri^.  Richardson  number  associated  with  temperature  differences. 
Equation  1 11-66 

Rig, Rig  Continuous  and  instantaneous  release  Richardson  numbers 
respectively 

Ri"  Richardson  number  associated  with  density  differences 

corrected  for  horizontal  turbulence  velocity 

Ri|  Richardson  number  associated  with  density  differences 

corrected  for  convective  scale  velocity 

Ri*  Richardson  number  associated  with  density  differences. 

Equation  III -56 

Sc  Schmidt  number 

Sh  Sherwood  number 

StH  Stanton  number  for  heat  transfer 

StM  Stanton  number  for  mass  transfer 

Sy  horizontal  concentration  scaling  parameter  (m) 

Sz  vertical  concentration  scaling  parameter  (m) 

Szq  Sz  at  the  downwind  edge  of  the  source  (x  =  1/2)  (m) 

2 

Sz0  value  of  SzQ  when  (tt  R  Q*)  is  a  maximum  (m) 
m 

T  temperature  associated  with  source  blanket  enthalpy  (K) 

T  .  temperature  associated  with  averaged  enthalpy  (K) 

C « L 

T$  surface  temperature  (K) 

TE  turbulent  kinetic  energy  (Nm) 

Tg  contaminant  storage  temperature  (K) 

t  time  (s) 
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characteristic  time  of  cloud  formation  (s) 
specified  time  (s) 

cloud  travel  time  to  the  maximum  distance  exposed  to  a 
given  concentration  (s) 

time  when  observer  i  encounters  downwind  edge  (s) 

time  when  observer  i  encounters  upwind  edge  (s) 

nondimensional  time  t/6 

ambient  average  velocity  (m/s) 

horizontal  or  frontal  entrainment  velocity  (m/s) 

effective  cloud  advection  velocity  (m/s). 

cloud  front  velocity  (m/s) 

velocity  of  observer  i  (m/s) 

local  velocity  in  j  direction  (m/s) 

average  transport  velocity  associated  with  (m/s) 

wind  velocity,  along  x-direction  (m/s) 

wind  velocity  measured  at  z  *  Zq  (m/s) 

horizontal  r.m.s.  turbulence  velocity  (m/s) 

internal  flow  out  of  gravity  current  head  (m/s) 

internal  flow  into  gravity  current  head  (m/s) 

wind  velocity  at  z  *  10  m  (m/s) 

friction  velocity  (m/s) 

characteristic  average  velocity  (m/s) 

cloud  volume  (m  ) 

heat  transfer  velocity  (0.0125  m/s)  in  Equation  1 1 1 -1 6 


3 

V.  initial  volume  (m  ) 

_  3 

specific  volume  of  air  (m  /kg) 

Vc  specific  volume  of  contaminant  (m  /kg) 

_  3 

Vw  uan-„  specific  volume  of  water  vapor  (m  /kg) 

w  j  Va pu r 

W  gas  plume  width  (m) 

Wq  rectangular  source  width  (m) 

w,  mass  fraction  of  air 

a 

w„  mass  fraction  of  contaminant 

c 

wg  vertical  entrainment  velocity  associated  with  HL  (m/s) 
w  total  entrainment  velocity  in  Equation  11-31  (m/s) 

w*  convective  scale  velocity  (m/s) 

Wg  entrainment  velocity  associated  with  (m/s) 

x. j (t)  x  position  of  observer  i  at  time  t  (m) 

xpl-  position  of  puff  center  due  to  observer  i  (m) 

x^  downwind  distance  where  gravity  spreading  terminates  (m) 

xy  virtual  point  source  distance  (m) 

xdn.  x  P°s‘it',on  downwind  edge  of  source  for  observer  i 

xup  x  P0S1^^0n  upwind  edge  of  source  for  observer  i 

x,y,z  Cartesian  coordinates  (m) 

Xq  downwind  edge  of  the  gas  source  (m) 

y  average  cloud  mole  fraction 

zR  surface  roughness  (m) 

Zg  reference  height  in  wind  velocity  profile  specification  (m) 


xxi  v 


•.•it k 


a  constant  in  power  law  wind  profile 

8  constant  in  oy  correlation  in  Equation  1 11-75 

8^  constant  in  Equation  11-37 

6'  coefficient  of  expansion  (K_1) 

0  .characteristic  time  scale,  vV®//gA.J 
r  gamma  function 

y  ratio  of  (p  -  pa)  /  cc 

Y-j  constant  in  Equation  III -75 

A  ratio  of  (p  -  p  )/p 

a 

AT  temperature  driving  force  (K)  (T  -  T  .  )  or  (T  -  T) 

S  C  )L  S 

A'  ratio  of  (p  -  pfl)  /  pfl 

A!  initial  ratio  of  (p  -  pa)/pa 

1  da 

6  constant  in  a  correlation  in  Equation  II 1-75 

6l  empirical  constant  (2.15)  in  Equation  III -31 
<5V  constant  (0.20)  in  Equation  11-59 
e  frontal  entrainment  coefficient  in  Equation  III-ll 

z,  collection  of  terms  defined  by  Equation  II 1-40 

X  Monin-Obukhov  length  (m) 

Xfus  latent  heat  of  fusion  for  water  in  Equation  A-4  (J/kg) 

Xvap  latent  heat  of  vaporization  for  water  in  Equation  A-4 

(J/kg) 

X'  modified  latent  heat  for  water  in  Equation  A-4  (J/kg) 
u  viscosity  (kg/m  s) 

q 

p  density  of  gas-air  mixture  (kg/m  ) 
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3 

ambient  density  (kg/m  ) 

3 

cloud  density  (kg/m  ) 

3 

initial  cloud  density  (kg/m  ) 

3 

vertically  averaged  layer  density  (kg/m  ) 

density  of  contaminant's  saturated  vapor  at  TQ  (kg/m  ) 

x-direction  dispersion  coefficient  (m) 

Pasqui 11 -Gifford  lateral  dispersion  coefficient  (m) 

z-direction  dispersion  coefficient  (m) 

2 

surface  shear  stress  (N/m  ) 

p 

shear  stress  at  the  cloud  top  (N/m  ) 
constant  of  Equation  I I -5 

function  describing  influence  of  stable  density 

stratification  on  vertical  diffusion,  Equation  III -54 
integrated  source  entrainment  function 
logarithmic  velocity  profile  correction  function 
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SUMMARY 


The  mathematical  modeling  methods  which  have  been  proposed 
for  prediction  of  dispersion  of  heavier-than-air  gases  are 
reviewed  and  critiqued.  The  phenomenology  of  atmospheric  dispersion 
of  heavy  gas  is  described,  and  a  simplified  representation  of  the 
general  heavy  gas  dispersion  scenario  is  proposed  which  involves 
three  distinct  flow  regimes: 

-Buoyancy-dominated  dispersion  (near  field) 

—Stably  stratified  shear  flow  (intermediate  field) 

—Passive  turbulent  dispersion  (far  field) 

Modeling  concepts  based  on  investigations  of  laboratory  gravity 
currents  are  used  to  describe  the  flow  and  dilution  processes  that 
characterize  the  buoyancy-dominated  dispersion  regime.  Measure¬ 
ments  of  the  spreading  and  dilution  of  heavy  gas  volumes  released 
suddenly  in  calm  (laboratory)  air  are  used  to  demonstrate  scaling 
methods  for  small  releases  from  35  to  530  liters,  and  the  laboratory 

O 

results  when  scaled  to  2000  m  volume  are  consistent  with  the 
experimental  data  from  the  Thorney  Island  2000  m3  Heavy  Gas  Trials. 

A  box  model  is  used  to  describe  the  laboratory  releases,  and  the 
gravity  spreading  and  air  entrainment  velocities  are  determined 
from  analysis  of  the  laboratory  experimental  data. 

Data  from  laboratory-stratified  shear  flow  mixing  experiments 
have  been  used  to  model  the  vertical  diffusion  of  heavy  gases  in 
the  atmospheric  constant  stress  layer.  The  modeling  concepts  used 
in  the  Shell  HEGADAS  (HEavy  GAs  Dispersion  from  Area  Sources)  model  have 
been  adapted  to  model  the  stably  stratified  shear  flow  and  passive 
turbulent  diffusion  regimes.  Several  important  modifications  to 
the  HEGADAS  model  have  been  made,  including  provision  for  heat 
transfer  and  convective  turbulence,  incorporation  of  additional 
laboratory-stratified  layer  mixing  data  in  the  vertical  entrainment 
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velocity  correlation,  and  improvements  in  the  method  used  for 
modeling  the  rate  of  heavy  gas  entrainment  (from  the  source)  by 
the  atmospheric  flow.  The  provision  for  air  entrainment  into  the 
gravity  spreading  gas  blanket  which  may  form  over  some  releases, 
as  well  as  heat  exchange  between  the  blanket  and  the  earth  surface, 
extend  the  applicability  of  the  model  to  a  much  wider  class  of 
heavy  gas  release  scenarios. 

The  interactive  computer  model  DEGADIS  (DEnse  GAs  Dispersion), 
which  is  proposed  for  incorporation  in  the  Coast  Guard  Hazard 
Assessment  Computer  System  (HACS),  has  been  used  to  simulate 
thirty-nine  experiments  which  represent  the  full  range  of  heavy  gas 
field  experiment  data  presently  available.  DEGADIS-predicted 
maximum  concentration  as  a  function  of  distance  is  compared  to  the 
maximum  reported  concentration  for  field  scale  releases  of  liquefied 
natural  gas  (LNG),  liquefied  petroleum  gas  (LPG),  and  Freon-12/air 
mixtures  from  the  Burro/Coyote,  Maplin  Sands,  and  Thorney  Island 
trials.  Based  upon  this  comparison,  the  variability  of  the 
distance  realized  to  a  concentration  level  of  5,  2-1/2,  or  1%  for  a 
given  release  is  quantified  based  on  the  predicted  distance.  The 
model  predicts  downwind  maximum  gas  concentration  decay  with 
distance,  to  the  hydrocarbon  flammable  range,  consistent  with  the 
full  range  of  experimental  data.  The  consistency  of  the  model 
predictions,  which  reflect  small  scale  laboratory  fluid  flow  and 
mixing  data,  with  the  results  of  the  wide  range  of  field  experimental 
releases  justifies  the  application  of  the  model  to  prediction  of 
dispersion  from  much  larger  releases. 

The  DEGADIS  model  can  readily  be  modified  to  incorporate 
better  mixing  and  heat  transfer  sub-models.  It  is  likely  that 
improvements  in  these  two  areas  can  be  expected  to  provide  the  most 
important  Information  for  improving  the  confidence  level  in  heavy 
gas  dispersion  prediction.  The  model  has  not  been  demonstrated  for 
application  to  the  prediction  of  dispersion  to  the  ppm  concentration 
range  which  characterizes  the  limits  of  toxic  gas  hazard.  However, 
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the  model  has  been  developed  to  be  consistent  with  the  existing 
data  base  on  passive  turbulent  dispersion  in  the  far  field,  and 
it  is  recommended  that  the  model  be  evaluated  against  experimental 
data  which  have  become  available  for  dispersion  of  ammonia  and 
nitrogen  tetroxide. 
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I.  INTRODUCTION 


Risk  of  accidental  release  of  heavier-than-air  gases  accompanies 
many  manufacturing,  storage,  and  transportation  operations.  Although 
increased  awareness  of  such  risks  has  been  fostered  by  debate  on  the 
risks  of  bulk  transport  of  large  quantities  of  flanmable  liquid 
gases  (General  Accounting  Office,  1978;  National  Research  Council, 

1980;  Office  of  Technology  Assessment,  1980),  which  in  general  are 
considered  to  be  heavier-than-air  gases  (HTAG),  many  other  chemicals 
handled  and  carried  in  bulk  can  produce  a  gas  or  aerosol  "cloud"  that 
is  more  dense  than  air  when  released  into  the  atmosphere  (Kaiser  and 
Walker,  1978;  Harris,  1978). 

Assessment  of  risk  attending  such  operations  involves,  in 
addition  to  estimation  of  the  probability  of  release,  prediction  of 
the  ensuing  atmospheric  dispersion,  since  such  dispersion  eventually 
results  in  dilution  of  the  gas  with  air  to  concentrations  which  are 
nonflammable  or  within  prescribed  toxicity  limits.  Consequently, 
a  prediction  of  the  location  of  the  "boundary"  of  such  clouds  (which 
may  be  defined  as  containing  gas  concentrations  above  a  prescribed 
lower  limit)  as  a  function  of  time  is  required  for  rational  risk 
assessment. 

Dispersion  of  HTAG  in  the  atmosphere  differs  importantly 
from  the  process  of  dispersion  of  trace  contaminants  in  the  atmosphere. 
The  theory  underlying  prediction  of  atmospheric  dispersion  of  trace 
contaminants  (pollutant  dispersion)  generally  assumes  that  the 
dispersion  is  the  result  of  the  turbulent  motion  that  characterizes 
the  atmospheric  boundary  layer.  The  presence  of  the  pollutant  is 
consequently  assumed  not  to  affect  the  atmospheric  flow  patterns; 
then,  the  problem  becomes  one  of  understanding  and  predicting  the 
atmospheric  boundary  layer  turbulence.  Although  the  characterization 
of  atmospheric  flow  suffers  from  the  general  limits  of  understanding 


turbulent  fluid  motion,  there  exists  a  fairly  well  developed 
theoretical  basis  for  prediction  of  the  dispersion  of  trace  contami¬ 
nants  in  the  atmosphere,  along  with  extensive  supporting  experimental 
data  derived  from  atmospheric  flow  measurements  (Pasquill,  1983). 

*  In  contrast,  the  release  of  large  quantities  of  HTAG  into  the 
atmosphere  may  alter  the  fluid  flow  pattern  in  the  atmosphere  in  the 
vicinity  of  the  release.  For  some  releases,  the  gravity-induced  flow 
and  the  resulting  interaction  with  the  atmospheric  flow  can  determine 
the  shape  and  extent  of  the  area  which  is  exposed  to  flammable  or 
toxic  gas  concentrations.  The  early  phases  of  cloud  formation, 
motion,  and  dispersion  following  release  of  HTAG  may  involve 
density-stratified  flows  which,  although  studied  widely  in  other 
contexts,  have  been  studied  from  the  perspective  of  atmospheric 
dispersion  of  HTAG  only  recently  (Turner,  1973). 

A  number  of  mathematical  modeling  techniques  for  predicting  HTAG 
dispersion  have  been  published.  The  models,  in  most  cases,  were 
developed  for  risk  assessment  studies  of  liquefied  gas  fuel 
importation  projects  in  the  U.S.  and  Europe.  When  applied  to  the 
prediction  of  the  maximum  downwind  distances  which  might  be  reached 
by  flammable  gas-air  mixtures  following  catastrophic  releases  of 
liquefied  natural  gas  from  ship  cargo  tanks,  order  of  magnitude 
differences  were  reported.  The  large  differences  in  predictions  of 
hazard  zones  extending  from  catastrophic  scenario  LNG  releases  were 
in  part  due  to  lack  of  standardization  of  release  and  atmospheric 
conditions  simulated,  but  also  reflected  different  modeling 
approaches  (Havens,  1977).  Since  1977,  extensive  laboratory  and 
field  experimental  programs  have  been  conducted  which  provide  data 
for  the  evaluation  of  HTAG  dispersion  models. 

This  study,  sponsored  by  the  U.S.  Coast  Guard  and  the  Gas 
Research  Institute,  attempts  to  provide  an  overall  assessment  of 
the  prediction  of  atmospheric  dispersion  of  HTAG  and  to  develop 
a  recommended  methodology  for  such  predictions  for  incorporation 
into  the  Coast  Guard  Chemical  Hazard  Response  Information 
System  (CHRIS). 


1.1  General  Description 


The  problem  is  prediction  of  the  dispersion,  by  mixing  with 
ambient  air,  of  a  mass  of  HTAG  released  during  a  finite  period  of 
time  in  the  atmospheric  boundary  layer.  The  initial  gas  cloud  may 
result  directly  from  (1)  an  aboveground  release  such  as  a  release 
of  pressurized  gases,  as  in  a  chlorine  tank  rupture,  or  (2)  indirectly 
from  a  ground-level  source  such  as  the  evaporation  of  a  released, 
liquid  as  in  the  case  of  a  spill  onto  the  sea  of  liquefied  petrol eum 
gas  (LPG)  or  natural  gas  (LNG).  The  cloud  may  be  initially  formed 
as  a  mixture  of  gas  or  gas-liquid  aerosol  and  humid  air.  Actual 
releases  will  always  be  of  finite  time  duration  and  will  have  a  finite, 
time-dependent  boundary  (extent). 

Appealing  to  generally  accepted  concepts  of  boundary  layer  flows, 
the  typical  development  and  movement  (Including  dispersion)  of  a 
heavy  gas  cloud  is  illustrated  in  Figure  1.1.  The  HTAG  source  is 
characterized  as  time-  and  position-dependent  with  the  source  mass 
rate  E  (x,y,t).  At  time  t,  a  HTAG  cloud  has  formed  over  the  source. 
Interaction  of  the  vertical  flux  of  gas  from  the  source  with  the 
atmospheric  flow  has  resulted  in  some  translation  downwind;  an  upwind 
movement  of  the  cloud  boundary  has  resulted  from  gravity-driven  lateral 
flow  of  the  heavy  cloud.  Characteristic  profiles  of  the  cloud 
boundary  at  times  t2  and  t3  after  completion  of  gas  injection  into 
the  atmosphere  are  depicted.  The  cloud  is  described  by  specification 
of  its  height  or  depth  H  (x,y,t),  local  concentration  c  (x,y,z,t), 
density  p  (x,y,z,t),  enthalpy  h  (x,y,z,t),  and  velocity  Uj(x,y,z,t). 
Characteristic  vertical  profiles  of  c,  H,  p,  and  u,  which  reflect  the 
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assumption  of  concentration  decrease  via  air  entrainment  through  the 
top  and  front  boundaries  of  the  cloud  (entrainment  velocities  wg  and 
ug  respectively)  as  well  as  shear  stress  at  the  cloud  top  and  bottom 
surfaces  (t^  and  tq  respectively),  are  illustrated. 

The  local  velocity  u^.  may  be  the  result  of  density-driving 
forces  which  cause  a  spreading  motion,  or  it  may  reflect  the 
atmospheric  flow  into  which  the  gas  is  injected.  In  actual  releases, 
some  combination  of  these  effects  will  be  important. 


The  HTA6  cloud  formation  illustrated  in  Figure  1-1  reflects  a 
cloud  formation  phase  whose  duration  is  small  compared  to  the  travel 
time  to  the  maximum  distance  exposed  to  the  concentration  of 
interest.  This  would  be  the  case,  for  example,  if  the  cloud  were 
formed  by  rapid  evaporation  from  a  rapidly  spilled  cryogenic  liquid. 

In  such  cases,  the  development  of  the  cloud  is  three-dimensional  and 
highly  transient,  and  a  general  model  of  the  momentum,  mass,  and 
energy  transfers  is  complex  and  difficult.  Several  mathematical  models, 
which  are  applications  of  different  numerical  computational 
procedures  for  solving  the  balance  equations  for  momentum,  energy, 
and  mass  with  initial  and  boundary  conditions  descriptive  of  the 
processes,  have  been  proposed  and  are  being  evaluated  (Havens,  1979, 
1983).  Such  models  are  complex  and  require  significant  computer 
solution  time  (and  cost).  Furthermore,  there  remain  important 
questions  regarding  the  validity  of  many  assumptions  made  in  such 
models  (i.e.  turbulence  closure  models)  which  require  further 
evaluation. 
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A  more  often  adopted  approach  to  modeling  the  dispersion  of  HTAG 
clouds  has  been  to  consider  the  limiting  cases  where  the  cloud 
formation  phase,  of  characteristic  duration  tf,  is  very  small  or  very 
large  with  respect  to  the  time  of  cloud  travel  to  the  maximum  distance 
exposed  to  the  concentration  of  interest,  tt-  When  t^  »  tt»  a 
stationary  "plume"  representation  of  the  HTAG  cloud  becomes  applicable. 
When  tf  <<  tt>  an  instantaneous  source  representation  is  indicated. 

In  either  case,  it  is  further  observed  that  a  typical  heavy 
gas  dispersion  scenario  involves  three  more  or  less  distinct  regimes 
of  fluid  flow.  Following  release,  especially  for  rapid  release  of 
a  large  quantity  of  heavy  gas,  a  cloud  having  similar  vertical  and 
horizontal  dimensions  (near  the  source)  may  form.  The  initial 
behavior  of  such  a  cloud  is  relatively  independent  of  the 
characteristics  of  the  ambient  wind  field  until  the  strength  of  the 
buoyancy-driven  flow  (slumping  and  lateral  spreading)  decreases 
sufficiently  that  the  cloud  motion  begins  to  be  determined  by  the 
ambient  atmospheric  flow.  When  the  cloud  notion  begins  to  be 
determined  by  the  atmospheric  flow,  the  dispersion  process  can  be 


described  as  a  stably  stratified  plume  (or  cloud)  embedded  in  the 
mean  wind  flow.  As  the  dispersion  proceeds,  the  stable 
stratification  due  to  the  heavy  gas  decreases  until  the  process 
can  be  represented  as  a  neutrally  buoyant  plume  (or  cloud)  in  a 
neutral  or  stratified  mean  wind  flow.  The  three  regimes, 

--  buoyancy-dominated  dispersion 
--  stably  stratified  shear  flow 
--  passive  dispersion  due  to  atmospheric  turbulence 
which  may  overlap  and  be  present  in  various  degrees  in  different 
heavy  gas  dispersion  scenarios,  must  be  accounted  for  if  a  model 
is  to  be  generally  applicable.  The  specific  treatment  of  each  of 
these  flow  regimes  in  earlier  models,  as  well  as  the  methodology 
used  to  provide  transition  between  the  regimes,  is  quite  varied 
and  explains  in  large  part  the  differences  observed  when  the 
various  models  have  been  applied  to  the  same  heavy  gas  dispersion 
scenario  (Havens,  1977,  1979). 

For  an^  atmospheric  dispersion  problem,  the  dispersion  of  the 
contaminant  is  passive  if  the  dispersion  is  entirely  the  result  of 
the  preexisting  atmospheric  turbulence,  and  the  contaminant  does  not 
affect  the  mean  flow  field  or  the  turbulence.  Consequently, 
conventional  air  pollutant  passive  dispersion  techniques  (Pasquill, 
1983)  become  applicable  to  the  description  of  the  latter  phases  of 
dilution  of  an  originally  HTAG  cloud.  A  method  for  delimiting  the 
zone  in  which  the  dispersion  process  is  non-passive,  i.e.  where  the 
flow  and  dispersion  processes  are  affected  by  the  presence  of  the 
HTAG,  is  needed  for  two  reasons.  First,  a  determination  of  whether 
the  release  need  be  described  as  a  HTAG  release  or,  more  simply,  as 
a  passive  pollutant  is  in  order.  Second,  experimental  data  must  be 
viewed  in  relation  to  the  relative  importance  of  the  gravity-driven 
flow  vs.  atmospheric  turbulence-dominated  processes  which  disperse 
the  gas. 


1.2  Buoyancy-Dominated  Dispersion 

Rapid  release  of  a  large  quantity  of  HTA6  may  result  in  an 
initial  cloud  having  similar  vertical  and  horizontal  dimensions.  The 
immediate  behavior  of  such  initially  "compact  clouds"  (Fay,  1984)  is 
controlled  by  the  ensuing  buoyancy-driven  flow  (slumping  and  lateral 
spreading).  Consequently,  the  initial  cloud  behavior  is  relatively 
independent  of  the  characteristics  of  the  ambient  wind  field  until 
the  strength  of  the  buoyancy-driven  flow  decreases  to  levels  at  which 
the  cloud  motion  begins  to  be  controlled  by  the  ambient  atmospheric 
flow. 

For  HTAG  releases  with  initially  similar  vertical  and  horizontal 
length  scales  there  is  now  conclusive  evidence  that  the  rapid 
gravity-driven  flow  which  ensues  results  in  large  scale  turbulent 
structures  which  effect  considerable  dilution  of  the  cloud  (Picknett, 
1978;  Hall  et  al.,  1982;  Meroney  and  Lohmeyer,  1982;  Havens  and 
Spicer,  1983;  Spicer  and  Havens,  1984).  Since  this  initial 
turbulent  motion  can  in  some  conditions  result  in  dilution  of  the 
cloud  by  a  factor  of  ten  to  one  hundred,  it  must  be  accounted  for 
in  predictions  of  HTAG  dispersion.  Furthermore,  an  understanding 
of  the  dilution  processes  present  during  tliis  phase  is  necessary 
to  accurately  predict  dispersion  of  large  scale  releases  in  low 
winds,  particularly  for  gas-air  mixture  flammability  levels,  which 
for  hydrocarbons  are  generally  of  order  one  percent. 

The  spreading  motion  that  follows  release  of  a  HTAG  volume  with 
similar  vertical  and  horizontal  length  scales  is  a  classical  "gravity 
current,"  characterized  by  the  gravity-driven  intrusion  of  one  fluid 
(the  HTAG)  into  i  less  dense  fluid  (the  atmosphere).  Gravity 
currents  are  formed  in  many  natural  situations  (Simpson,  1982), 
including  thunderstorm  outflows,  sea  breeze  fronts,  and  cold  front 
movements  in  the  atmosphere,  and  a  variety  of  ocean  currents  driven 
by  temperature,  salinity,  or  suspended  solids-induced  density 
differences. 


A  comprehensive  review  of  laboratory  gravity  current  studies  and 
their  relation  to  gravity  currents  in  the  atmosphere  and  ocean  has 
been  recently  published  by  Simpson  (1982).  Gravity  currents  formed 
by  instantaneous  release  of  a  fixed  volume  of  denser  fluid  under  a 
layer  of  less  dense  fluid  pass  through  three  flow  regimes.  There  is 
first  a  "slumping  phase,"  during  which  the  current  is  retarded  by  the 
counterflow  of  the  bulk  fluid  into  which  the  intrusion  occurs;  this 
slumping  phase  prevails  until  the  ratio  of  the  depth  of  the  density 
current  to  the  bulk  fluid  is  about  0.075  (Huppert  and  Simpson,  1980). 
For  fixed  volume  releases  in  the  atmosphere,  this  phase  is  extremely 
short  and  can  be  neglected  from  our  consideration.  In  the  immediately 
ensuing  flow  phase  the  dynamics  of  the  current  are  governed  by 
inertial  and  buoyancy  forces.  In  this  phase,  the  buoyancy  and 
inertial  forces  may  be  roughly  considered  to  be  in  balance,  and  the 
motion  of  the  bulk  fluid  is  unimportant.  Finally,  there  follows  a 
viscous  phase,  wherein  the  buoyancy  force  is  roughly  balanced  by 
viscous  forces.  The  viscous  phase  of  gravity  currents  resulting 
from  the  fixed  volume  release  of  a  HTAG  is  characterized  by  the 
virtual  absence  of  further  dilution  by  turbulent  mixing; 
accordingly,  we  concentrate  on  the  buoyancy-inertial  phase  of 
such  flows. 

During  the  inertial -buoyancy  phase,  the  motion  of  the  current 
can  be  considered  quasi -stationary,  and  the  distinguishing  features 
of  the  flow  near  the  gravity  current  head  are  represented  in  Figure 
1.2. 


Figure  1.2.  The  head  of  a  steady  gravity  current  (Simpson 
and  Britter,  1979,  and  van  Ulden,  1983). 


Characteristic  features  of  the  flow  process  illustrated  in  Figure 
1.2  are  (van  Ulden,  1983): 

1.  At  the  leading  edge  a  "head"  is  present  with  a  depth  H-| 
about  twice  that  of  the  current  behind  the  head,  Hh. 

2.  An  elevated  stagnation  point  is  present,  below  which  an 
insignificant  flux  Q-j  of  ambient  fluid  is  entrained. 

3.  Behind  the  head  a  wake  region  is  present  in  which 
significant  mixing  occurs.  An  associated  significant 
internal  flow  in  the  head  is  present  with  a  denser  fluid 
layer  flowing  into  the  head  from  behind  with  a  velocity 

=  0.2  times  the  velocity  of  the  front,  and  a  mixed 
layer  flowing  away  from  the  head. 

4.  The  front  velocity  relative  to  the  flow  being  intruded 
into  is  well  described  by  the  relation 

Uf  ■  CE(9A'Hh)1/2  (1-1) 

with  CE  =  1.15  +  0.05,  based  on  the  data  from  Schmidt 
(1911),  Benjamin  (1968),  Fannelop  et  al.  (1980),  and 
Huppert  and  Simpson  (1980). 

The  use  of  the  relation  given  in  Equation  1-1  to  describe  the 
lateral  gravity-driven  spreading  of  a  fixed  volume  release  of  HTAG, 
coupled  with  the  representation  of  the  gas  cloud  as  a  vertically 
oriented  cylinder  whose  radius  and  height  change  as  a  result  of 
gravity  spreading  and  air  entrainment  across  the  outer  surfaces 
of  the  cloud  has  formed  the  basis  for  several  HTAG  dispersion 
models  (van  Ulden,  1984;  Germeles  and  Drake,  1975;  Fryer  and 
Kaiser,  1979;  Cox  and  Carpenter,  1979;  Fay,  1980;  Picknett,  1981; 
Eidsvik,  1980;  Fay  and  Ranck,  1981;  and  Fay,  1984). 

These  models  have  been  reviewed  by  Havens  (1982)  and  Webber 
(1983).  Meroney  and  Lohmeyer  (1982)  and  van  Ulden  (1979,  1983) 
have  proposed  bulk  models  for  describing  gravity  currents  which 
account  for  the  acceleration  phase  of  the  current  which  precedes 
the  quasi -steady  inertial -buoyancy-controlled  phase  depicted  in 
Figure  1.2.  A  summary  of  these  modeling  techniques  drawn  upon  in 
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formulating  a  model  for  the  buoyancy-controlled  phase  is  oven  in 
Section  II.  An  extensive  laboratory  investigation  of  the  gravity 
spreading  and  dilution  of  right  circular  cylindrical  volumes  of 
HTAG  released  in  calm  air,  from  which  is  derived  pertinent  data  for 
modeling  the  gravity  spreading  and  air  entrainment  during  the 
buoyancy-controlled  phase  of  such  releases,  is  described  and 
discussed  in  Volume  II  of  this  report. 

1.3  Stably  Stratified  Shear  Flow 

An  intermediate  phase  of  the  dispersion  process  is  cnaracterized 
by  its  similarity  to  a  wide  variety  of  naturally  occurring  flow 
processes  in  which  a  stably  stratified  plume  is  embedded  in  a  mean 
flow.  Such  plumes  are  expected  to  differ  importantly  from  neutrally 
buoyant  (passive)  plumes  in  three  significant  ways. 

1.  A  lateral  (crosswind)  gravity-driven  mean  flow  will  persist 
until  the  negative  buoyancy  of  the  cloud  has  been  reduced 
(by  entrainment)  to  render  the  mean  flow  velocity  small 
compared  to  an  appropriate  lateral  turbulent  velocity 
characteristic  of  the  atmospheric  flow. 

2.  The  density  interfaces  which  define  the  plume  boundaries 
should  act  to  damp  turbulent  mixing  and  consequently  reduce 
vertical  mixing  (air  entrainment)  into  the  HTAG  cloud. 

3.  Because  such  plumes  are  characterized  by  a  small  vertical 

to  longitudinal  ratio,  dilution  occurs  primarily  as  a  result 
of  vertical  mixing. 

The  ground  level  lateral  extent  of  a  negatively  buoyant  plume 
embedded  in  a  turbulent  boundary  layer  has  been  investigated  by 
Britter  (1980).  Considering  the  steady  release  of  a  gas  with  density 

3 

c,  at  a  volumetric  rate  Q  m  /s,  from  a  circular  source  of  diameter  D 
at  ground  level  into  a  wind  field  with  average  velocity  ua,  the  cross- 
section  of  the  cloud  at  the  downwind  edge  of  the  source  is 
approximated  as  rectangular  with  width  D  and  height 

H  =  Q  /  Wnu  .  If  it  is  assumed  that  the  lateral  gravity-driven 

U  a 


m 
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spreading  is  described  by  Equation  1-1  and  the  gas-air  mixture 
buoyancy  is  preserved,  the  cloud  width  as  a  function  of  downwind 
distance  is  then  given  by 


r  ' 

1/2 

3  CE 

9  A'  Q 

x  +  W2/2 

3 

ua  j 

(1-2) 


and  the  rate  of  lateral  spread  with  respect  to  downwind  distance  is 


dx 
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(1-3) 


The  lateral  gravity-driven  spreading  rate  is  indicated  to  be  a 
function  of  the  length  scale  gA'  Q/u^,  hence  the  lateral  spreading 

Q 

rate  in  a  given  wind  field  is  determined  by  the  buoyancy  flux  QgA'. 
The  lateral  extent  of  the  cloud  is  then  dependent  in  the  near-field 

3 

on  the  two  length  scales  D  and  g  A'  Q/u_. 

u 

Britter  has  experimented  with  the  injection  of  a  dense  fluid 
(saline  solution)  through  the  floor  (area  source)  of  a  water  flume. 
The  flume  was  0.76  m  wide,  10  m  long,  and  the  mean  flow  in  the  flume 
was  characterized  as  fully  turbulent  with  Reynolds  number,  based  on 
bulk  velocity  and  flume  depth,  greater  than  5,000  and  a  friction 
velocity  equal  to  the  bulk  velocity  divided  by  20.  The  dense  fluid 
source,  which  was  flush  with  the  flume  floor,  was  a  sintered  bronze 
distributor  plate  with  diameter  5  cm,  placed  in  the  centerline  of  the 
flume.  The  ranges  of  variables  studied  were 


0  <  Q  <60  cm3/s 
0  <  u,  <40  cm/s 

a 

0  <  g  A'  <  200  cm/s^ 


VAs 
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Britter's  experiments  indicated  that  the  buoyancy-driven  lateral 

velocities  were  greater  than  the  friction  velocity  of  the  mean  flow 
3  4 

for  xua  /  Q g A  »  1 .3  x  10  .  It  follows  that  the  transition  to 

a 

passive  behavior  could  be  expected  to  be  complete  when 
4  Q  Q  A 1 

x  >  1 . 3  x  1 0  z-y -  and  the  pi ume  shoul d  be  passive  from  the  source 
ua 

uj  D 

when  - >  1000. 

Q  g  A' 

Fluid  mixing  across  the  density  interface  which  characterizes 

a  dense  plume  embedded  in  a  turbulent  boundary  layer  was  reviewed 

by  Turner  (1973).  There  is  general  agreement  that  entrainment 

across  such  plume  boundaries  is  correlated  with  a  bulk  Richardson 

2  2 

number  for  the  plume  defined  as  Ria  =  g  A'  H/u.  or  Ri*  =  gA'H/u*. 

The  turbulence  which  effects  mixing  in  a  plume  flowing  over  a 
solid  surface  and  under  a  turbulent  flow  layer  may  originate  in  the 
shear  stress  at  the  bottom  solid  surface  or  at  the  density  interface. 
For  thin  layers  "intruding"  under  a  layer,  the  bottom  stress  is 
probably  most  important,  and  this  case  probably  best  represents 
the  behavior  of  HTAG  plumes  in  the  atmosphere. 

McQuaid  (1976),  Kantha  et  al .  (1977),  and  Lofquist  (1960)  have 
reported  experimental  studies  from  which  entrainment  velocity 
correlations  can  be  derived  for  bulk  Richardson  numbers  encompassing 
the  range  of  interest  for  HTAG  dispersion  in  the  atmosphere.  McQuaid 
reported  measurements  of  the  vertical  dispersion  of  carbon  dioxide 
(c  '  o,  =  1.52)  introduced  as  a  floor-level  line  source  in  a  fully 

Q 

developed  rectangular  channel  flow  in  a  wind  tunnel  of  width  0.3  m, 
height  0.9  m,  and  working  length  5.5  m.  The  mean  velocity,  u  ,  in 

a 

the  tunnel  (=  0.8  umav)  ranged  from  0.82  to  3.52  m/s  with  a 
friction  velocity  (in  the  absence  of  any  HTAG)  reported  as 
u*  =  0.0412  u_.  The  channel  Reynolds’  number  (based  on  the  mean 

a 

velocity  and  the  channel  hydraulic  diameter)  ranged  correspondingly 
4  5 

from  2.5  x  10  to  1.1  x  10  .  Defining  an  average  plume  depth  via  the 


relation  H 


- - - ,  where  urt-c  is  the  effective  advection 

cmax  UEFF 


velocity,  the  entrainment  velocity  of  air  into  the  two-dimensional 
plume  becomes  w^  =  u£pp d Hppp /  dx.  Britter  analyzed  McQuaid's 
data  and  reported  the  correlation  for  entrainment  velocity  with 
Ri*  shown  in  Figure  1.3. 

Kantha  et  al .  measured  turbulent  entrainment  at  the  interface 
in  a  stably  stratified  two-layer  fluid  with  constant  surface  stress 
applied  at  the  free  (top)  surface.  The  rate  of  increase  of  the 
mixed  layer  depth  with  time  was  measured  photographically  using 
a  dye  initially  distributed  in  the  top  (lighter)  layer.  The 
Richardson  number  (Ri*  *  Hj_  g  a'/  u*)  is  constant  in  such  an  experiment. 
One  value  of  the  entrainment  velocity,  defined  as 


s 

| 

$ 


we  = 


(1-4) 


was  obtained  for  a  given  value  of  u*  and  A  (i.e.  for  a  given 
Richardson  number).  Kantha  et  al.  reported  measurements  of  wg  for 
30  <  Ri*  <  1000.  In  order  to  compare  the  entrainment  velocities 
obtained  by  Kantha  et  al . ,  defined  by  Equation  1-4  with  the  data 
from  McQuaid,  must  be  related  to  the  mixed  layer  depth  H^. 
Assuming  that  the  entrained  material  is  represented  as  being 
distributed  vertically  with  a  Gaussian  concentration  profile  and 
that  H^pp  is  defined  as  the  height  where  the  concentration  is 
one-tenth  that  of  the  maximum  value,  HL  =  2.15  H^pp.  The  data  of 
Kantha  et  al .  expressed  as  w^  =  dH^pp/dt  is  shown  in  Figure  1.3. 

Lofquist  (1960)  also  measured  turbulent  entrainment  across  a 
stable  density  interface  formed  between  a  turbulent  shear  layer  of 
salt  water  and  a  quiescent  pure  water  layer  above  and  presented  the 
entrainment  velocity  as  a  function  of  a  Froude  number 

Fr2  «  - —  (1-5) 


9  A'  H 


where  Hr  is  the  hydraulic  radius  of  the  rectangular  channel  in  which 
the  experiments  were  conducted.  The  data  presented  by  Lofquist  have 
been  re-expressed  in  the  form 


Correlation  of  entrainment  velocity  with  bulk  Richardson  number 
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d  HPjrp 

w;»^-.f(R1.)  (1-6) 

and  included  in  Figure  1.3. 

The  combined  data  of  McQuaid,  Kantha  et  al.,  and  Lofquist 
shown  in  Figure  1.3  were  curve-fitted  to  give 

w^/u*  =  k(l  +  a)/ [0.88  +  0.099  Ril*04 

+  1.4  x  10‘25  Ri*-7]  (1-7) 

where  a  is  the  power  for  a  power  law  wind  velocity  profile 
(Equation  I I I -50) .  This  functionality  is  obtained  if  the  vertical 
turbulent  diffusivity  is  given  by 


Kz  *  k  u*z  /  <j>  (Ri*) 


and  the  vertical  concentration  distribution  is  of  the  form 


c(x,z)  *  cc(x)  exp 


* 

1+a 

z 

S  (X) 

l  2  J 

as  developed  in  Section  III.  For  curve-fitting,  a  value  of  0.2 
was  used  for  a.  The  limiting  value  (Ri*  *  0)  of  w^  /  u*  *  l  /0.88 
incorporates  the  ratio  of  eddy  diffusivities  for  heat  (thermal 
energy)  and  momentum,  Km/KH,  and  k  =  0.35,  as  determined  ty 
Businger  et  al .  (1971)  from  measurements  in  the  neutral  atmospheric 
surface  stress  layer.  The  eddy  diffusivity  for  mass  K£  is  assumed 
equal  to  that  for  heat  (energy)  transfer  under  the  conditions 
reported  by  Businger  et  al. 

The  effect  of  vertical  density  stratification  on  vert' cal 
mixing  shown  by  the  data  is  consistent  with  a  passive  limit 
for  Wg  /  u*  of  0.4  suggested  by  the  solution  of  the  diffusion 
equation  (Equation  1-19)  with  a  constant  wind  velocity  and  mass 
diffusivity  given  by 


1.4  Passive  Dispersion  due  to  Atmospheric  Turbulence 

At  some  distance  from  the  source,  gas  will  be  sufficiently 
dilute  to  justify  its  consideration  as  a  trace  material,  and  the 
dispersion  can  be  modeled  using  approaches  developed  for  atmospheric 
pollutant  dispersion  (Hanna  et  al.,  1982;  Pasquill,  1983). 

Considering  quasi -steady  releases  at  ground  level  in  the 
atmospheric  flow  field  from  a  Lagrangian  point  of  view,  the  crosswind 
(lateral)  and  vertical  spread  of  a  contaminant  occurs  as  a  result  of 
small  scale-distorting  processes.  However,  the  trajectories  of 
successive  sections  of  the  plume  may  be  irregularly  displaced  by 
larger  scale  fluctuations  in  the  flow  (plume  meander).  Thus  the 
time  average  concentration  experienced  downwind  of  a  point  source 
not  only  diminishes  with  distance  from  the  source  but  also  varies 
with  time.  This  property  of  the  time-mean  concentration  results 
from  dispersive  motions  with  length  scales  larger  than  the  plume 
cross-section.  It  follows  that  the  time-mean  crosswind  distribution 
from  a  source  may  not  be  described  adequately  as  a  diffusion  process 
with  the  gradient  transfer  hypothesis  because  of  the  displacement 
of  the  plume  by  such  larger  scale  fluctuations  in  the  flow.  However, 
it  is  assumed  that  vertical  spread  from  a  ground  level  source  is 
controlled  by  dispersive  motions  whose  scale  is  restricted  by  the 
Dresence  of  the  boundary  (and  therefore  are  small  with  respect  to  the 
vertical  spread),  and  vertical  dispersion  from  ground  level  sources 
has  been  widely  modeled  using  the  gradient  transfer  hypothesis  in  the 
diffusion  equation.  Vertical  dispersion  in  the  surface  layer  from 
ground  level  sources  is  conventionally  modeled  using  the  gradient 
transfer  hypothesis  combined  with  application  of  similarity  principles 
developed  by  Monin  (1959)  and  Batchelor  (1964)  and  extended  for 
stratified  flow  by  Gifford  (1962). 

The  vertical  velocity  gradient  responsible  for  maintaining  the 
stress  in  the  surface  layer  is  assumed  to  be  a  function  of  the 
surface  stress  Tq,  the  density  p,  and  the  height  z.  From 
dimensional  principles 
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E  u*  /  kz 


(1-10) 


Equation  1-10  defines  the  friction  velocity  u*,  and  the  experimentally 
determined  value  of  k  (von  Karman  constant)  is  variously  reported  as 
0.35  to  0.40.  Integration  of  Equation  1-10  gives  the  logarithmic 
velocity  profile 


In  -L 


(1-11) 


where  zQ  is  an  integration  constant  which  represents  the  height  at 
which  the  velocity  is  taken  to  be  zero. 

Equation  I -1 1  is  valid  for  neutrally  stable  flow.  If  there  is 
an  upward  or  downward  flux  of  heat  (qQ)  through  the  surface  layer, 
vertical  density  variation  will  result,  and  the  density  stratification 
results  in  modification  of  the  flow  field.  The  vertical  momentum  flux 
is  then  the  result  of  buoyancy  flux  as  well  as  mechanical  shear;  the 
flux  of  buoyancy  is 


r  - 

F°-?V 


(1-12) 


where  8'  is  the  expansion  coefficient  of  the  fluid  in  the  flow  field. 
A  length  scale  which  quantifies  the  relative  importance  of  shear- 
related  and  buoyancy-related  vertical  momentum  transfer  is  defined 
by  combining  the  buoyancy  and  momentum  fluxes: 


A  *  u;/F0 


(1-13) 


where  X  is  the  Monin-Obukhov  length.  For  qg  negative  (heat  and 
momentum  transfer  in  the  same  direction),  X  is  positive,  and  the 
associated  stable  stratification  of  the  fluid  layer  tends  to  suppress 
turbulence.  If  qg  is  positive,  an  unstable  density  stratification 
results  (X  <  0)  ,  tending  to  increase  vertical  turbulent  transfer. 

A  neutrally  stable  fluid  layer  corresponds  to  infinite  X. 
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From  dimensional  arguments,  a  modified  form  of  Equation  I -1 0 
is  derived  to  describe  the  vertical  velocity  profile  in  a  stratified 
shear  layer 


dux  _  u* 
dz  ”  k?  N 


(1-14) 


with  <t>M  determined  from  experimental  measurements  of  vertical 


momentum  transfer.  A  widely  used  analytic  representation  of  <j> 


M 


based  on  field  wind  profile  measurements  and  proposed  by  Businger 
et  al.  (1971)  is: 


1  - 16  f 


-1/4 


X  <0  (unstable) 


(1-15) 


X  >  0  (stable) 


(1-16) 


The  vertical  eddy  diffusivity,  defined  as  the  ratio  of  momentum 
flux  to  the  mean  velocity  gradient, 


Tn  I  P 


du. 


dz 


(1-17) 


becomes,  using  Equation  1-14  and  noting  that  tq  *  p  u* 
ku*z 


Si 


M 


(1-18) 


For  two-dimensional,  steady  state,  vertical  diffusion  in  the 
surface  stress  layer,  the  diffusion  equation  is 


uf£-i 
x  3x  3z 


v  9c 
c  3Z 


(1-19) 


where  «c  is  the  mass  eddy  diffusivity.  It  is  common  practice  to 
invoke  Reynolds'  analogy  between  mass  and  momentum  transfer 
processes  to  assume  Kc  =  K^,  although  there  is  experimental 
evidence  that  the  eddy  diffusivities  for  mass  and  heat  or  thermal 
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energy  (K^)  are  greater  than  KM.  Businger  et  al .  (1971)  showed 
that  the  ratio  was  about  1.3  in  neutrally  stable  flow 

(A  -  °°) ,  decreases  to  near  1.0  with  increasing  stability 
(z/A  -+  2.0)  and  increases  to  about  2.5  for  decreasing  stability 
(z/A  -  -2.5). 

The  application  of  Equation  1-19,  with  suitable  descriptions 
of  ux(z)  and  Kc (z ) ,  for  modeling  vertical  dispersion  in  the  surface 
layer  of  the  atmosphere  has  been  demonstrated  (Hanna  et  al.,  1982; 
Pasquill,  1983). 
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II.  REVIEW  OF  AVAILABLE  MODELS 

A  number  of  mathematical  techniques  for  predicting  HTAG 
dispersion  have  been  published.  These  models  can  be  (somewhat 
arbitrarily)  classified  in  three  categories. 

1.  Box  models  which  represent  the  initial  development 
of  a  cloud  (in  the  case  of  an  instantaneous  release) 
or  a  cross-sectional  slice  of  a  cloud  (in  the  case  of 
a  steady  continuous  release)  as  a  uniformly  mixed 
volume.  The  shape,  thermodynamic  properties,  and 
position  of  the  cloud  are  modeled  using  correlations 
derived  for  the  velocity  of  density  intrusions  and 
fluid  entrainment  across  density  interfaces.  This 
model  type  often  incorporates  a  transition,  usually 
to  a  Gaussian  model,  to  describe  passive  dispersion 
(controlled  by  atmospheric  turbulence)  of  the  gas  in 
the  far  field. 

2.  K-theory  models  which  assume  constitutive  relations 
between  turbulent  fluxes  and  the  gradients  in  mean 
variables  velocity,  temperature,  and  concentration, 
coupled  with  the  time-averaged  incompressible  flow 
balance  equations  for  mass,  momentum,  and  energy  to 
predict  the  temporal  and  spatial  variation  of  cloud 
properties. 

3.  There  are  several  models  which  involve  greater 
simplification  of  the  equations  of  motion,  energy,  and 
mass  than  is  found  in  the  K-theory  models.  Some  still 
require  solution  of  partial  differential  equations  to 
predict  cloud  state  variables  such  as  the  models  proposed 
by  Zeman  (1980),  Rosenzweig  (1980),  and  Fannelop  (1980) 
and  refinements  thereof,  such  as  the  SLAB  model  published 
by  Lawrence  Livermore  National  Laboratories  (Morgan  et 
al.,  1983a)  which  are  one-dimensional  time-dependent  models 
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designed  to  describe  diffusion  and  gravity-driven  flow 
of  a  HTAG  release  in  the  atmosphere.  Several  models 
(te  Riele,  1977;  Flothmann  and  Nikodem,  1980;  Colenbrander, 
1980)  assume  power  law  concentration  profiles  of  the 
form 


c(x,y,z)  =  cc(x)  exp 

and  some  couple  this  assumption  with  a  K-theory  represen¬ 
tation  of  turbulent  mass  transfer  within  the  cloud 
(Colenbrander,  1980). 

The  K-theory  and  SLAB  models  were  not  considered  suitable 
for  incorporation  in  an  interactive  computer  system  such  as  HACS 
and  will  not  be  described  here.  It  is  appropriate  to  note  that  the 
SLAB  model  and  FEM3,  a  K-theory  model  published  by  Lawrence  Liver¬ 
more  National  Laboratories  (Chan  et  al.,  1982),  utilize  concepts 
to  describe  the  effect  of  stable  density  stratification  on  vertical 
dispersion  which  are  similar  to  those  used  in  other  models  such  as 
Colenbrander  (1980),  Eidsvik  (1980).  and  the  model  recomnended 
herein. 

Vertical  entrainment  of  air  into  a  gas  layer  is  expected  to  be 
a  function  of  the  stabilizing  effect  of  the  density  gradient 
relative  to  the  shear  flow  characterized  by  vertical  turbulence. 

If  the  HTAG  flow  is  viewed  as  being  superimposed  on  the  local 
atmospheric  flow  and  convection-induced  turbulence  is  neglected, 
a  characteri stic  vertical  turbulence  velocity  can  be  represented  as 
a  friction  velocity  of  the  flow,  where 

u*  =  (1  cp)1/2u  (II -2 ) 

Cp  is  a  surface-friction  drag  coefficient  and  u  is  a  characteristic 
wind  velocity  ua,  or,  as  represented  in  some  models,  the  weighted 
vector  sum  of  the  gravity  spreading  velocity  and  the  characteristic 
wind  velocity 


f  ' 
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u2  =  (c*uf)2  +  u2  ( 1 1-3) 

The  vertical  density  stratification  of  the  flow  is  measured  by  a 
form  o f  the  overall  Richardson  Number 


where  w  is  given  by 

2  2  2 

w  =  u*  +  nw*  ( 1 1  -  5 ) 

w*  represents  the  convection  scale  velocity,  and  n  is  an  empirical 
constant  (Zeman  and  Tennekes,  1977). 

Classical  boundary  layer  analysis  (Pasquill,  1983;  Turner,  1973) 
suggests  that  in  the  constant  stress  layer  of  the  atmosphere,  the 
vertical  entrainment  velocity  should  be  proportional  to  the  friction 
velocity  in  the  absence  of  stratification  and  inversely  proportional 
to  the  Richardson  number  for  stratified  flow. 

wg  -  c-|  u*  for  Ri  -*>  0  (1 1-6) 

we  =  c2  U*/Rl  for  Ri  >:>  1  ( 1 1 -7) 

These  equations  have  been  demonstrated  for  boundary  layers  adjacent 
to  the  earth's  surface. 

Velocity  shear  is  also  produced  at  the  interface  between  a 
spreading  HTAG  layer  and  the  overlying  atmospheric  flow. 

Turbulent  mixing  across  such  a  surface  can  result  from  breaking 
waves  associated  with  instabilities  generated  at  the  interface 
(Turner,  1973).  Some  HTAG  cloud  models  have  provided  for  air 
entrainment  via  this  mechanism;  the  entrainment  is  assumed 
proportional  to  the  difference  in  velocity  across  the  cloud  top. 

For  a  spreading  "top  hat  profile"  gas  layer,  such  entrainment  is 
represented  by 


where  u  is  the  average  cloud  top  velocity.  In  the  following,  some 
approaches  to  air  entrainment  modeling  are  summarized. 

II.l  Box  Models 

A  number  of  models  have  been  proposed  in  this  category;  however 
it  is  possible  to  consider  general  forms  which  have  been  used  to 
model  gravity  spreading  and  air  entrainment  into  the  homogeneous 
volume  or  cloud.  This  model  type  is  further  categorized  based  on 
the  methodology  used  to  describe  the  gravity-driven  HTAG  flow.  In 
earlier  models,  the  dynamics  of  the  HTAG  flow  was  simply  described 
as  a  quasi -steady  density  intrusion,  reflecting  the  assumption  of 
exchange  of  cloud  potential  energy  for  lateral  kinetic  energy  of  the 
flow;  this  category  is  referred  to  as  a  hydrostatic  approximation 
box  model.  Refinement  to  include  a  momentum  balance  on  the  cloud 
spreading  process  which  describes  more  realistically  the  acceleration 
from  rest  as  well  as  the  ensuing  slumping  and  lateral  spreading 
process  have  appeared  in  some  models. 

Since  these  models  assume  uniform  thermodynamic  properties, 
their  use  is  restricted  to  the  representation  of  two  types  of  HTAG 
clouds,  illustrated  in  Figure  II.l. 

1.  An  instantaneously  released  quantity  of  gas  with  a  known 
initial  volume,  usually  represented  as  a  vertically 
oriented  cylinder  whose  radius  and  height  change  as  a 
result  of  gravity  spreading  and  air  entrainment  across 
the  outer  surfaces  of  the  cloud.  The  cylindrical  cloud 
is  usually  assumed  to  be  translated  with  wind. 

2.  A  steady  plume  with  rectangular  cross-section  (axis  along 
the  wind  direction)  whose  thermodynamic  properties  vary 
with  downwind  position  and  whose  width  and  height  change 
as  a  result  of  gravity  spreading  (assumed  lateral  only) 
and  air  entrainment  across  the  outer  surfaces. 


Wind  Direction 


Source 
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Cloud  shape  for  Instantaneous  Release 
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Cloud  shape  for  Continuous  Releases 


Figure  II. 1  Box  model  cloud  representation 
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In  either  case,  the  models  require  analytical  expressions  for 
the  spreading  velocity  (i.e.  the  velocity  of  the  cloud  edge)  and 
the  entrainment  of  air  at  the  cloud  boundaries. 

1 1. 1.1  Hydrostatic  Approximation  Box  Models 

All  proponents  of  this  approach  have  modeled  the  velocity  of 
the  cloud  edge  using  a  density  intrusion  relation 


uf  *  CE 


->1/2 


P  “  P. 


P  J 


=  C£[gAH] 


1/2 


(II-9) 


or 
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Cc[gi'H]’/2 


(11-10) 


Equation  1 1 -10  follows  from  Equation  1 1 -9  when  p  =  pfl.  The  density 
of  the  cloud  is  affected  by  energy  transfer  from  the  cloud 
surroundings  as  well  as  the  entrainment  of  air  into  the  cloud.  No 
consideration  of  the  different  approaches  to  modeling  the  surface- 
to-cloud  heat  transfer  has  been  included  here.  It  is  probable, 
however,  that  such  heat  transfer  may  importantly  affect  the  dispersion 
of  cryogenic  gases  such  as  LNG;  some  investigators  have  indicated 
that  air  entrainment  may  be  significantly  enhanced  due  to  convection¬ 
generated  turbulence  (Eidsvik,  1980;  Fay  and  Ranck,  1981).  Both  of 
these  effects  are  discussed  in  Section  IV. 

Although  box  model  treatments  of  air  entrainment  have  given 
widely  disparate  results,  most  of  the  differences  are  attributable 
to  differences  in  interpretation  of  the  sparse  data  base  used  to 
determine  empirical  entrainment  coefficients.  Box  models  can  be 
represented  as  requiring  entrainment  velocities  which  are  applied  in 
the  general  form 
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Vpa  =  weAT  +  ue  AF 


(11-11) 


where  wgAj  and  ug  Ap  are  vertical  and  horizontal  entrainment  rates 


represented  as  the  product  of  a  characteristic  area  and  velocity. 
Because  the  HTA6  clouds  are  generally  shallow  and  grow  laterally 
due  to  gravity  forces,  the  second  term  in  Equation  I I -11  is  normally 
assumed  to  be  important  only  in  the  gravity  spreading-dominated  stages 
of  the  cloud  development.  The  horizontal  entrainment  velocity  has 
most  often  been  modeled  by 


ue  =  e  uf 


(II-12) 


or 


e  u. 


ue  =  uf(t=0T 


(11-13) 


Equation  11-13,  proposed  by  Eidsvik  (1980),  reduces  to  Equation  11-12 
in  the  limit  t  -*■  0,  and  gives  an  effective  horizontal  entrainment 
velocity  which  decreases  to  zero  for  A-p/Ap  »  1.  Details  of  specific 
models  are  summarized  chronologically  in  the  following. 


van  Ulden  (1974)  characterized  the  HTAG  flow  with  the  Richardson 


number  Ri*  =  gAH  /  u*.  He  suggested  that  cloud  frontal  entrainment  is 


very  small  and  that  cloud  top  entrainment  during  the  gravity-dominated 

2 

spread,  which  he  defined  as  Ri*»  4/Cp,  is  negligible.  Thus 
van  Ulden's  model  of  the  gravity  spread  phase  indicates  minimal 


cloud  volume  growth  (and  dilution)  given  by 


V°.  1  $  ‘  ueAF  S  °-°5i.fAF 


van  Ulden  suggested  that  atmospheric  turbulence-induced  mixing  will 


begin  to  control  when  the  spreading  velocity  has  been  reduced  to  2  u* 

o 

(Ri*  =  4  /  C| )  and  thereafter  models  the  cloud  dilution-using  a 
classical  Gaussian  passive  dispersion  model.  It  appears  that 


van  Ulden's  (1974)  model  entrainment  parameters  were  based  primarily 
on  analysis  of  the  growth  of  a  heavy  Freon-air  cloud  following 
essentially  instantaneous  vaporization  of  1000  kg  liquid  Freon 
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dumped  into  water.  It  is  important  to  note  that  measurements  of 
the  initially  formed  visible  cloud  indicated  a  rapid,  initial, 
approximately  ten-fold  dilution  of  the  Freon,  and  van  Ulden  modeled 
the  subsequent  dispersion  only. 

Germeles  and  Drake  (1975)  neglected  cloud  frontal  entrainment 

and  modeled  top  entrainment  using  Equation  I 1-8.  Germeles  and  Drake 

suggested  c3  *  0.1  based  on  their  Froude  number  extrapolation  of 

Lofquist's  (I960)  data  for  entrainment  across  a  density/shear  interface 

between  a  salt  water  flow  and  overlying  quiescent  pure  water.  The 

flow  in  Lofquist's  experiments  was  three-dimensional  (the  flow  was  in 

a  channel  with  comparable  width  and  depth),  and  although  Lofquist's 

correlation  of  the  data  utilized  a  Froude  number  which  incorporated 

the  hydraulic  radius  of  the  heavy  layer,  this  appears  not  to  have  been 

accounted  for  in  Germeles'  and  Drake's  extrapolation.  Analysis 

indicates  an  order  of  magnitude  uncertainty  in  the  value  of  the 

entrainment  coefficient  of  0.1  inferred  from  the  Germeles-Drake 

extrapolation  if  the  geometry  of  the  flow  is  taken  into  account. 

Also,  the  top  entrainment  coefficient  of  0.1  inferred  by  Germeles 

and  Drake  is  sensitive  to  the  value  of  C£  in  the  front  spreading 

velocity  equation  (Equation  11-10).  If  Cc  =  1.0  is  assumed,  in 

1/2  t 

contrast  to  C^  =  2  as  suggested  by  Germeles-Drake,  their  extrapola¬ 
tion  gives  03  *  0.01  at  the  front  velocity.  Finally,  the  extrapolated 
value  of  Cg  =  0.1  would  appear,  in  any  case,  to  be  applicable  only 
to  the  top  area  near  the  cloud  front  (since  it  reflects  the  Froude 
number  at  the  front)  and  should  be  smaller  for  those  areas  of  the 
cloud  top  which  are  spreading  at  lower  velocity. 

Cox  and  Carpenter  (1980)  recommended  modeling  frontal  entrainment 
using  Equation  11-12  with  e  =  0.6  and  top  entrainment  using  the 
following  expressions 

w  =  0.15  u,  for  Ri"  -  0 

6  '  (11-15) 

W g  =  0.1  u-j  /  Ri"  for  Ri"  »  1 


where  u-j  is  the  horizontal  r.m.s.  turbulence  velocity  of  the  wind  flow 
and  Ri"  is  defined  as 


R1 "  • 


(11-16) 


u^  is  assumed  to  be  proportional  to  u*,  so  that 


u,  =  —  u* 
1  u*  * 


(11-17) 


Values  of  u-j  /  u*,  attributed  to  Pasquill  (1974),  of  3.0,  2.4,  and 

1.6  for  very  unstable,  neutral,  and  very  stable  atmospheric  flows 

respectively,  were  used.  The  turbulence  scale  length  L  is  correlated 

with  height  above  ground  (cloud  depth)  and  atmospheric  stability 

using  data  from  Taylor  et  al.  (1970).  Using  the  suggested  values  for 

u-|  /  u*  along  with  estimates  of  L  /  h  from  Taylor  et  al . ,  Equation 

11-15  can  be  rewritten  in  the  form  of  Equations  I I -6  and  II -7, 

along  with  a  stability-dependent  proportionality  constant,  as 

follows:  we 

Vertical  Entrainment 
_ Velocity 

Pasquill  Stability  Class  Ri"-*-  0  Ri1' »  1 


0.45  u„ 
0.36  Uj 
0.24  u. 


3.0  u*/Ri" 
2.0  u*/Ri" 
0.7  u*/Ri" 


The  Cox  and  Carpenter  model  changes  to  a  Gaussian  profile  model  when 

dav 

the  cloud  spreading  velocity  decreases  to  u^  =  ua  where  oy(x) 
is  the  crosswind  dispersion  coefficient  for  point  source,  passive 
plumes  proposed  by  Pasquill  and  Gifford  (Pasquill,  1974). 

Fryer  and  Kaiser  (1979)  used  the  same  top  and  frontal 
entrainment  velocity  models  as  Cox  and  Carpenter  but  suggested 
transition  to  a  Gaussian  dispersion  model  when  either  of  the  following 
two  conditions  apply: 
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u^  <  2.14  ua  2  and  Ri  <  0.1 


or 


(11-18) 


A'  <  8.16  x  10 


-4 


Fay  (1980)  suggested  considering  both  wg  and  ug  proportional  to 
u^  during  the  gravity  flow-dominated  regime  (u^  »  u*),  so  that 


/  pa  c4ufAT  +  c5UfAF 


(H-19) 


Based  on  his  analysis  of  Ellison  and  Turner's  (1959)  laboratory 
experiments  with  steady  stratified  fluid  layers  flowing  above  or  below 
a  quiescent  fluid.  Fay  suggested  an  upper  limit  of  0.01  for  c^,  and 
based  on  his  analysis  of  Simpson  and  Britter's  (1979)  gravity 
intrusion  experiments,  he  suggested  a  value  for  c5  of  0.05. 

Picknett  (1981)  correlated  the  concentration  and  visible  cloud 
dimension  data  from  a  series  of  instantaneous  releases  on  land  of 
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approximately  40  m  of  Freon-air  mixtures  ranging  in  initial  relative 
density  (p/p,)  from  about  2.0  to  5.0.  Picknett  fitted  the  large 

a 

initial  entrainment  inferred  from  the  data  using  the  frontal  entrain¬ 
ment  relation  (Equation  11-12)  with  e  =  0.82.  He  modeled  top 
entrainment  using  Equation  1 1 -7  with  C2  =  0.15.  Picknett  terminated 
the  gravity  spreading  phase,  and  transitioned  to  a  Gaussian  passive 
dispersion  model  when  the  Richardson  number  of  the  spreading  layer 
decreased  to  Ri*  =  7. 

Eidsvik's  (1980)  top  hat  model  does  not  transition  to  a  Gaussian 
passive  dispersion  model.  Instead,  the  model  uses  air  top  entrainment 
relations  incorporating  Equation  I I -7  (dominant  during  the  early 
gravity  spread  phase  when  Ri*>>  1)  and  Equation  1 1 -6  (dominant  during 
the  later  phase  of  the  cloud  dilution  when  density-driven  flows  have 
subsided,  as  indicated  by  Ri*-*-  0).  Incorporating  Eidsvik's  suggested 

1 

model  parameter  values  and  using  a  different  Ri*,  his  equation  for 
vertical  (top)  entrainment  is 
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we  = 


0.39  w 

1  +  Ri!/  19.7 


(11-20) 


which  in  the  absence  of  thermal  convection  effects  (Eidsvik  included 
recommended  measures  of  convective  turbulence  velocities  which  are 
not  considered  at  this  point)  is  given  by  Equation  I I -5  with  n  *  0.29. 

The  friction  velocity  u*  is  determined  using  Equation  I 1-2  and 
a  surface  drag  coefficient  along  with  a  characteristic  velocity  u 
calculated  using  Equation  1 1 -3;  the  constant  c*  is  dependent  on  the 
geometry  of  flow  (c*  *  2/3  for  axi symmetric  and  1/2  for  one-dimensional 
spreading).  The  Richardson  number  Ri*  is  then  defined  as 


(11-21) 


Using  Equations  II-2,  3,  and  9,  the  Richardson  number  Ri*'  can  be 
expressed  as 

..2 


£  (1/2  cf)((c*uf)Z  +  uZ) 


In  the  limit  as  /  ug  <<1, 


(11-22) 


=  0.39  w 


(11-23) 


uf  «  ufl 


which  is  consistent  with  the  assumed  dependence  of  w£  on  u*  given  by 
Equation  I I -6 . 

For  Uf  >>  ug,  however,  wg  approaches  (for  c^  <<1.0,  =  1.3, 

and  c*  *  2/3) 


1.35  cj/2  Uf 


uf  »  ua 


(11-24) 
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Eidsvik  cites  recommended  values  for  c^  of  2  x  10  for  flow  over 
water  and  1.4  x  10"^  for  flow  over  land.  For  this  range,  and  with 
=  1.3,  c*  =  2/3,  Wg  given  by  Equation  11-24  varies  between 
0.0001  u.f  and  0.002  uf.  It  should  be  noted  that  Eidsvik's  limiting 
(as  uf  »  u  )  entrainment  expression  of  Equation  11-24,  with  a 

'  a  _4 

coefficient  of  order  10  to  10  ,  infers  a  different  entrainment 

mechanism  from  that  of  Germeles'  and  Drake's  suggested  entrainment 
relation  (Equation  I I -8) .  Eidsvik's  model  neglects  entrainment  due 
to  shear  at  the  interface  (he  indicates  that  Kelvin-Helmhol tz 
instability  might  provide  a  mechanism  for  such  entrainment,  but 
assumes  it  to  be  effectively  damped  out  because  of  the  strong  density 
stratification  across  the  interface)  whereas  Germeles'  and  Drake's 
suggested  entrainment  is  associated  with  the  shear  at  the  interface. 

In  a  subsequent  review  of  the  box  model  formulation  for  HTAG 
dispersion.  Fay  and  Ranck  (1981)  correlated  the  experimental  data 
from  HTAG  wind  tunnel  releases  by  Hall  (1979)  and  Meroney  (1982) 
and  field  experiments  by  Feldbauer  et  al .  (1972),  van  Ulden  (1974), 
Picknett  (1978),  and  Koopman  et  al.  (1979,  1981)  with  a  zero  frontal 
entrainment  (cf  =  0)  and  top  entrainment  velocity  given  by 
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C1  c2 


Rif 


+  c 
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(11-25) 


Fay  and  Ranck  reported  correlation  of  their  model  with  the  isothermal 
HTAG  release  data  they  analyzed  which  indicated  c^  *  2.5  and  C£  =  0.5. 
They  stated  that  the  model  correlation  with  data  for  nonisothermal 
gases,  which  was  primarily  LNG  experimental  release  data,  indicated 
c-|  =  2.5  and  c2  =  5.0. 

*  Finally,  Fay  (1984)  has  suggested  that  for  initially  compact 

/Clouds,  i.e.  where  vertical  and  longitudinal  initial  dimensions  are 
-jjiifci.'lar,  there  is  rapid  mixing  during  the  period  following  release 
during,  which  the  buoyancy-dominated  flow  is  controlling.  Based 
•primarily  on  his  analysis  of  the  no-wind  HTAG  release  data  of  Hall 
(1982),  Meroney  and  Lohmeyer  (1982),  and  Havens  and  Spicer  (1983), 


Fay  proposed  modeling  the  entrainment  during  the  buoyancy-dominated 
flow  period  with  the  relation 


(11-26) 


which  reflects  the  assumption  of  a  constant,  global  entrainment  rate 
with  a  constant  cloud  height.  For  releases  in  wind.  Fay  proposes 

Vea  =  0-6  v5/6  (9A')1/2+c6(ttR2)  u*/Ri*  (11-27) 

with  Cg  =  0.6,  as  representative  of  the  dispersion  in  the  previously 
mentioned  no-wind  releases  as  well  as  near  instantaneous  wind  tunnel 
releases  reported  by  Hall  (1979,  1982)  and  Meroney  and  Lohmeyer 
(1982). 
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II. 1.2  Non-Hydrostatic  Approximation  Box  Models 


van  Ul den  (1979),  Meroney  and  Lohmeyer  (1982),  and  van  Ulden 
(1983)  have  proposed  methods  for  modeling  the  gravity  slumping/ 
spreading  phase  of  instantaneously  released  HTA6  clouds.  These 
methods  incorporate  a  momentum  balance  to  describe  the  acceleration 
from  rest  and  subsequent  motion  of  a  HTAG  flow. 

van  Ulden  (1979)  analyzed  a  nonstationary  gravity  current 
resulting  from  an  instantaneous  release  of  a  mass  M  of  HTAG. 
Described  here  is  van  Ulden's  analysis,  applied  to  the  axisymmetric 
lateral  spreading  of  an  instantaneously  released  "cloud." 

An  energy  balance  on  a  control  volume  sufficiently  large  to 
contain  all  influences  of  the  spreading  cloud  requires  that  the 
total  energy  remain  unchanged  during  the  spreading  process: 


d  PE  .  d  KE  .  d  TE  .  d  IE 

~ST  +  “Tr“  + 


dt 


dt 


dt 


(11-28) 


In  Equation  11-28,  PE  and  KE  represent  the  potential  and  kinetic 
energy  of  the  cloud,  respectively.  TE  represents  the  turbulent 
kinetic  energy,  and  IE  represents  the  irreversible  conversion  of 
kinetic  (and  turbulent)  energy  to  internal  thermal  energy. 

If  the  volume,  mass,  and  density  of  the  cloud  are  initially 
V.j ,  M.. ,  and  p.. ,  respectively,  and  the  density  of  surrounding  air  is 


Pa»  the  buoyancy  of  the  gas  is  g 
energy  is 


pi  -  pa 


Pa 


V. ,  and  the  potential 


PE  =  (M.  -  pa  V.)  g  H/2  =  g  A M H/2 


(11-29) 


Since  AM,  the  amount  by  which  the  initial  gas  mass  exceeds  that  of 
an  equivalent  volume  of  air,  is  constant,  it  follows  that 


d  PE 
~3T 


1  nAMdH 

7  gAMdt 


(11-30) 


If  it  is  assumed  that  the  spreading  gas  cloud  assumes  the  shape  of 
a  right  cylinder. 
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(11-31) 


The  first  term  on  the  R.H.S.  of  Equation  11-31  represents  the 
downward  velocity  of  the  cylinder  top  surface  which  would  occur  if 
the  volume  of  the  cylinder  remained  constant,  and  the  second  term 
represents  the  upward  velocity  of  the  top  surface  due  to  an  increase 
in  volume  by  entrainment  of  air. 

The  kinetic  energy  of  the  cylindrically  shaped  gas-air  mixture 
volume  is  given  by 


KE 


K 


r  H  (u(r) ;dr 


(11-32) 


Continuity  requirements  dictate  that  u(r)  *  ^uf  »  50  that  Equation 
11-32  gives 


KE(t)  *  {  Mu2 


(11-33) 


and  it  follows  that 


dKE  .  1  M  „  duf  .  1  .  2  dM 
■ar  ‘  2Muf  “dT  +  4  uf  dt 


(II-34) 


Since  dM/dt  =  p  Q  *  ttR2p  w  ,  Equation  11-34  can  be  written  as 

d  ut  a  i 


dKE 

“3F 
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1  u  uuf  .  1  2  M  Dawt 

2  Muf  ~5t+  4  ufM  TTT 


(11-35) 


A  balance  equation  for  energy  associated  with  turbulence  in  the 
flow  is  proposed: 

=  turbulent  energy  produced  by  shear,  SP  (11-36) 

-  turbulent  energy  transformation  into  potential 
energy  due  to  entrainment,  "buoyant 
destruction,"  BD 

-  turbulent  energy  dissipation  into  thermal 

d  IE 

energy, 

It  is  assumed  that  the  production  of  fluid  turbulence  by  shear  occurs 
mainly  at  the  leading  edge,  and  contributions  from  friction  at  the 
cloud  lower  surface  and  upper  surface  behind  the  leading  edge  are 
neglected.  Since  the  advancing  cloud  front  encounters  undisturbed 
air  at  a  rate  proportional  to  2ttRHu^,  and  the  undisturbed  air  must 
be  accelerated  to  a  velocity  of  the  order  of  uf,  van  111  den  suggests 
that  the  rate  of  turbulent  energy  production  by  shear  be  represented 
as 

SP  «  B1(2TrRHuf)(l  pa  UJ)  =  y  £  uf  (II-37) 

Assuming  that  the  buoyant  destruction  term  in  Equation  11-36 
is  just  equal  to  the  potential  energy  gain  by  entrainment  of  air  and 
substituting  Equations  11-35,  36,  37  into  Equation  11-28,  the  energy 
balance  equation  for  the  "cloud"  is  written  to  obtain  an  equation  for 
the  front  velocity 


duf  T  Da  wt 
~St  +  2uf  pTT 


2qApH  ^a  ,  2 

p  R  R  p  f 


=  0 


(11-38) 


If  air  entrainment  is  neglected  and  the  motion  is  assumed  stationary 

du.f 

(  "fo  =  0),  Equation  11-38  reduces  to 
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(fI-39) 


and  the  hydrostatic  spreading  velocity  model  is  retrieved. 

Meroney  and  Lohmeyer  (1982)  used  Equation  11-38  as  developed 
by  van  Ulden,  modified  to  reflect  the  assumption  that  the  turbulent 
energy  production  rate  at  the  cloud  edge  scales  with  the  cloud  density 
rather  than  air  density: 


i  .e. 


3  3 

uf  pa  uf 

SP  =  $i  M-p-  instead  of  B-|  M 


(11-40) 


Meroney  and  Lohmeyer  assumed  that  during  the  gravity  spreading  phase 
of  the  HTAG  motion,  entrainment  would  occur  at  the  top  and  front 
surfaces  of  the  cloud  as  in  Equation  I I -11  and  further  represented 
wg  and  ug  as  proportional  to  the  front  velocity  uf 


ue  =  C1  uf 


(11-41) 


we  =  c2  uf 


Then 


=  2c^ttRHu^  +  c^ttR^u^ 


(11-42) 


Replacing  the  total  entrainment  rate  (Equation  11-31),  Equation  11-38 
was  expressed  by  Meroney  and  Lohmeyer  as 


du 


dt 


f  +  vv 
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c2  .  C1 
2  H  R 


(11-43) 


Meroney  and  Lohmeyer  conducted  a  series  of  instantaneous  releases  of 
hemi cylindrical  volumes  of  Freon-air  mixtures  in  calm  air  and,  based 
on  the  analyses  thereof,  suggested  the  following  values  for  the 

l  l 

constants  in  Equation  11-43:  B-j  =  0.9,  c^  =  0.1,  C2  =  0.1. 


van  Ulden  (1983)  proposed  a  new  bulk  model  which  also  incorporates 
a  balance  equation  for  momentum  to  describe  the  development  of  an 
entraining  gravity  current  in  rectangular  coordinates.  The  current 
is  divided  into  a  "head"  and  "tail"  with  characteri sti c  lengths  and 
velocities  as  illustrated  in  Figure  II. 2.  Considering  a  radial 
ax i symmetrical  current,  the  front  moves  with  velocity  uf(t)  at 
position  R(t) .  Tne  tail  region,  which  extends  to  r  =  R^,  advances 
with  the  layer-averaged  velocity  u^(t).  The  height  (or  depth)  of 
the  tail  and  head  regions  are  and  respectively.  In  the  following 
van  Ulden's  development  is  adapted  for  a  cylindrical  coordinate  system. 
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Figure  II. 2  The  unsteady  gravity  current  (van  Ulden,  1983) 
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The  radial  component  momentum  balance  on  the  gravity  spreading 
cloud  is  represented  as 

dP  _  P  /static  pressure  force,  due  to  negative 
dt  p  [buoyancy  of  the  gravity  current  / 

+  Fp  [dynamic  force,  equal  to  the  sum  of  the 

I  drag  force  on  the  head  of  the  current  and 
]the  lift  force  arising  from  asyiranetry  of 
(the  ambient  flow  around  the  current  head 

+  F  (force  due  to  acceleration/deceleration  '• 
[reaction  by  the  ambient  fluid  / 

(11-44) 

The  static  pressure  force,  obtained  by  integrating  the  static 
pressure  over  the  boundary  of  the  current,  is 

Fp  =  (jgApHt)(2rRHt)  =  tt  g  A  p  R  H2  (11-45) 

Neglecting  the  shear  stress  at  the  bottom,  the  dynamic  force  on 
the  current  is  the  sum  of  the  drag  force  on  the  head  of  the  current 
and  the  lift  force  that  arises  due  to  asymmetry  in  the  ambient  flow 
around  the  head.  The  drag  force  is  represented  as 

FD  =  *  T  pa  uf  (27rRhavHh^  =  "  avdvirRHhpa  uf  (H-46) 
where  dv  is  an  effective  drag  coefficient  and  the  constant  ay  is  an 
empirical  ratio  of  the  average  head  depth  H-j  to  Hh  (ay  =  H-j  /  Hh ) . 

The  horizontal  acceleration  reaction  Fa  is  approximated  by  the 

a 

reaction  to  an  accelerating  elliptical  cylinder  with  an  aspect  ratio 
H/R  (Batchelor,  1967): 

Fa„  =  "  <JT  (kl  (II-47) 

K 

and  the  vertical  acceleration  reaction  is  represented  as 


1 


it  (k2pa*RH  uf> 


(11-48) 


where  k-j  and  kg  are  coefficients  of  order  one. 
constant.  Equations  11-47  and  11-48  give 


Using  a  single 


Fa  =  'v"0. 


d(R  IT  uf) 


(11-49) 
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Using  Equations  11-45,  46,  and  49,  the  momentum  balance  (Equation 
11-44)  becomes 


dP  2  2  d(R  H  uJ 

3t  *  -  avdv”caRHhuf  -  e»  71  pa  - — - 

dt 


(11-50) 


As  before  (1979),  van  Ulden  assumed  that  the  potential  energy 
decrease  due  to  slumping  of  the  cloud  is  offset  by  the  production 
of  kinetic  energy,  which  through  the  action  of  shear,  is  partly 
transformed  to  turbulent  kinetic  energy.  Part  of  the  turbulent 
kinetic  energy  is  transformed  back  into  potential  energy  due  to 
entrainment  of  air  by  the  cloud,  van  Ulden  assumes  that  this 
"buoyant  destruction"  of  kinetic  energy  is  proportional  to  the  rate 
of  production  of  turbulent  kinetic  energy,  and  appealed  to  the 
analysis  of  Simpson  and  Britter  (1979)  to  argue  that  the  turbulent 

3 

kinetic  energy  production  rate  should  scale  as  ttp  HRu^.  Then, 


}gApH  =  £7rpaHRu^ 


(II-51) 


which  can  be  written 
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(11-52) 


where  e  is  an  empirically  determined  coefficient. 


The  volume  integral 
R 


V  = 


h(r,t)rdr 


(11-53) 


0 

where  h(r,t)  is  to  be  expressed  in  terms  of  Hh  and  H-t,  and  the 
momentum  integral 


■/ 


P  =  2ir  pu(r,t)  h(r,t)  rdr  * 


pt  *  ph 


(11-54) 


are  then  approximated  with  separate  analyses  of  the  head  and  tail  of 
the  current. 

In  the  tail  of  the  current,  the  shallow  water  equations  are 
assumed  applicable.  It  is  assumed  that  the  shape  of  the  current 
is  quasi -stationary  in  time,  and  the  layer-averaged  density 
difference  is  assumed  horizontally  uniform.  It  follows  that  the 
volume  and  momentum  of  the  tail  are  given  by 
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Vt  “  *Rh  (Ht  +  V  /2*  (II-55) 

Pt  *  1°  *lHt  +Hh)  ”Rh  T  (II-56> 

A  momentum  balance  for  the  head  region,  assuming  quasi -steady 
state,  indicates  that  the  static  and  dynamic  pressure  forces  on 
the  head  should  be  balanced  by  the  net  flux  of  momentum  due  to 
flow  into  and  out  of  the  head.  The  static  pressure  and  drag  are, 
respectively 

Fp  *  (l/2sAoHh)(2iTRhHh)  * 

Fd  ■  -dv(l/2Fau2f)[2»Rh<\Hh>] 

2 

=  "  2.A,  P,Ux  7T  R.  H. 

v  v  a  t  n  n 

Near  the  surface,  the  inward  flow  (u^  in  Figure  1.2)  carries 

momentum  into  the  head,  while  the  return  flow  (u,  in  Figure  1.2) 

0  1 

carries  momentum  out  of  the  head.  Assuming  u3  =  u^,  *  £Hh, 

and  u4  =  <5vuf>  the  momentum  flux  into  the  head  is  approximately 

Qh  *  6vpauf  C27rRhHh]  (II'59) 


(11-57) 


(11-58) 


with  *  0.2,  dy  =  0.64;  Equation  11-60  then  specifies  the  head 
velocity  boundary  condition.  The  volume  of  the  head  is  determined 
by  assuming  that  the  head  length  scales  with  H-j ,  and  it  follows 
that 


R  -  Rh  =  byH1  (11-61) 

where  bv  is  an  empirical  constant.  The  volume  of  the  head 
becomes 

\  -  *\bv  (R  +  Rh)Hh  (H-62) 

If  the  layer-averaged  velocity  is  assumed  to  increase  linearly 
with  r,  it  follows  that 


Along  with  the  definition  of  u^. 


(11-63) 


(11-64) 


(11-65) 


Equations  11-50,  52,  60,  61,  63,  and  65  can  be  solved  to  determine 
o,  Ht,  Hh,  V,  Ph,  and  Pt> 
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I I. 1.3  Comparison  of  Box  Model  Predictions  for  No-Wind  Releases 

Figure  II. 3  illustrates  the  substantial  differences  in 
dilution  of  a  1  m  cylinder  (H/D  =  1 )  of  Freon-12  gas  (sp.  gr. 

=  4.19)  instantaneously  released  in  calm  air,  as  predicted  by  the 
models  proposed  by  van  Ulden  (1974),  Germeles-Drake  (1975), 
Cox-Carpenter  (1979),  Fryer  and  Kaiser  (1979),  Eidsvik  (1980), 

Fay  (1980),  Meroney  and  Lohmeyer  (1982),  van  Ulden  (1983),  and  Fay 
(1983)  are  described  in  the  previous  two  sections. 

Figure  II. 3  also  summarizes  the  cloud  average  concentration 
vs.  cloud  radial  extent  obtained  by  spatial  integration  of  the 
radial  and  vertical  cloud  concentrations  determined  from  the 
laboratory  calm-air  heavy  gas  releases  conducted  under  this 
contract  (Volume  II,  this  report).  Finally,  the  average  cloud 
concentration  decay  with  distance  predicted  with  a  model  adapted 
from  van  Ulden  (1983),  which  incorporates  empirical  constant 
values  to  fit  the  experimentally  determined  cloud  average 
concentration  decay  with  distance  (or  time),  is  also  shown  in 
Figure  1 1. 3.  This  model,  which  has  been  incorporated  in  the 
DEGADIS  model  recommended  herein,  will  be  described  in  detail  in 
Section  III. 


II. 2  Advanced  Similarity  Models 

Experimental  concentration  data  from  passive  ground  level 
source  plume  dispersion  experiments  are  generally  well  described  as 
Gaussian  (Hanna,  1982;  Pasquill,  1983).  The  spatial  distribution 
of  concentration  from  a  point  source  is  represented  as 


where  cc(x)  is  the  centerline  concentration  at  ground  level  (z  =  0) 
and  Oy,  cz  are  the  standard  deviations  of  the  concentration 
distribution  in  the  y  and  z  directions.  Correlations  of  ay  and  oz 
as  functions  of  distance  downwind  of  the  plume  source  and  atmospheric 
stability  have  been  developed  based  on  a  combination  of  theoretical 
and  field  experimental  data.  The  most  frequently  used  correlations 
of  dispersion  data  in  terms  of  o  and  az  are  based  on  field 
atmospheric  tracer  dispersion  experiments  conducted  during  the  1950s 
and  1960s.  The  Project  Prairie  Grass  experiments  (Haugen,  1959) 
involved  concentration  measurements  at  distances  to  about  1  km 
downwind  of  near  ground  level,  point  source,  steady  releases  of  a 
tracer  material;  they  are  the  basis  for  the  oy  and  c>2  correlations 
developed  by  Pasquill  (1961)  and  adapted  by  Gifford  (1961,  1968, 

1976)  into  the  form  given  in  Figure  II. 4.  Analytical  formulas  for 
ay  and  a2,  based  on  the  data  from  the  Prairie  Grass  experiments  as 
well  as  experimental  data  from  passive,  elevated  (z  ■  108  m)  tracer 
release  experiments  conducted  at  Brookhaven  National  Laboratory, 
have  been  proposed  by  Smith  (1968),  Briggs  (1973),  and  Pasquill 
(1975).  Laboratory  data  from  stably  stratified  shear  flow  experiments 
also  support  a  Gaussian  description  of  vertical  concentration  profiles 
(McQuaid  (1976)). 

The  assumption  of  a  Gaussian  similarity  form  for  vertical  and 
horizontal  concentration  profiles  in  a  HTAG  plume  was  proposed  by 
te  Reile  (1977)  and  further  developed  by  Colenbrander  (1980,  1983). 
Both  considered  a  model  for  atmospheric  dispersion  of  HTAG  emitted 
from  an  idealized  rectangular  source  of  width  2B  and  length  L  as 
illustrated  in  Figure  I I. 5.  The  gas  source  center  is  at  x  =  y  *  z  *  0 
with  x  representing  downwind  direction,  and  y  and  z  horizontal 
(crosswind)  and  vertical  directions  respectively.  Power  law  profiles 
for  the  velocity  and  gas  concentration  are  assumed: 


Figure  II. 4.  Field  data  correlations  for  a  and  a  for  passive  dispersion 
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Figure  I I. 5. 


Form  of  assumed  concentration  profiles-- 
te  Reile  (1977)  and  Colenbrander  (1980). 


u •y,v*.»Vi 


:.v 


49 


The  concentration  profile  has  a  center  section  of  width  2b  which  is 
represented  as  dispersing  only  in  the  vertical  direction  to  account 
for  the  relatively  uniform  concentration  field  which  develops  over 
a  uniform  area  source.  Horizontal  diffusion  processes  are 
associated  with  gas  concentration  gradients  at  the  "edge"  of  the 
uniform  center  section.  As  a  result  of  the  similarity  assumptions 
in  Equations  11-67  and  11-68,  the  variables  cc»  b,  S^,  and  S2  are 
functions  of  downwind  distance  (x)  only  and  are  therefore  constrained 
by  ordinary  differential  equations. 

Both  te  Reile  and  Colenbrander  suggested  modeling  the  vertical 
and  horizontal  dispersion  scales  Sz  and  by  requiring  that  the 
concentration  profiles  satisfy  the  corresponding  two-dimensional 
diffusion  equations 
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u  Z—  m  JL. 
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with  turbulent  diffusivities  described  by 
k  u*  z 

Kz  *  7WJ 
and 


Ky  K0uxW  1 


(11-69) 


(11-70) 


(II-71) 

(11-72) 


Colenbrander  proposed  the  function  <p 

<t>(Ri*)  =  0.74  +  0.25  Ri*0,7  +  1.2  x  10”7  Ri*3  (11-73) 

based  on  his  analysis  of  HTAG  entrainment  velocities  measured  by 
McQuaid  (1976)  and  Kantha  et  al  (1977)  (see  Section  I). 
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Equation  11-70  models  the  horizontal  turbulent  diffusion  as  a 
gradient  transfer  process.  As  previously  stated,  such  a  representation 
is  questionable  because  the  horizontal  scales  of  turbulence  may  not  be 
small  compared  to  the  time  average  plume  width.  Consequently,  the 
determination  of  in  Equation  11-70  as  a  function  of  x  is  forced 
to  reflect  results  of  experimental  data  on  horizontal  dispersion 
of  passive  plurres  from  point  sources,  such  as  the  oy(x)  correlations 
developed  by  Pasquill  and  Gifford,  te  Reile  and  Colenbrander  assumed 
that  a  horizontal  turbulent  diffusivity  is  proportional  to  the  mean 
wind  velocity  and  to  the  "width"  of  the  cloud  raised  to  a  power  y 
where  the  width  is  defined  (see  Figure  1 1. 5)  to  be 


(11-74) 


Noting  that  S  =  /T  a  and  representing  the  a  vs.  x  correlation  of 

J  J  J 

Figure  II. 4  by  ay  =  5xe  where  6  and  3  are  constants  dependent  on 
the  Pasquill -Gifford  stability  class.  Equations  11-70,  11-72,  and 
11-74  determine  an  equation  for  Sy.  Colenbrander  modeled  the  lateral 
growth  of  the  plume  due  to  gravity  spreading  using  the  gravity 
intrusion  velocity  specification  of  Equation  I I -9: 


g(p ( cc )  -  Pa(z  *  o)) 

p(cc) 


(11-75) 
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with  the  effective  depth  H^pp  defined  by 
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(11-76) 
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Defining  a  mean  downwind  advection  rate  u^pp 
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(11-77) 


it  follows  that  the  downwind  lateral  growth  of  the  cloud  can  be 
modeled  by 


dw  =  1  dW 

dx  ~  Uppp  dt 


(11-78) 


Equations  which  determine  Syl  Sz,  and  b  as  a  function  of  x  are 
integrated  until  b  =  0  (denoted  by  x  *  xt).  At  that  time,  lateral 
gravitational  spreading  of  the  cloud  is  terminated,  and  Sy  at 
greater  distances  is  specified  by  the  relation: 


Sy-/T  o y  -  /T  «(x  +  xv)B  (11-79) 

where  xy  is  the  virtual  point  source  location  determined  from  the 
value  of  Sy  at  x  =  xt< 

Finally,  the  total  mass  flux  of  gas  for  a  steady  release  is 
independent  of  x  and  given  by 


2  cc(x)u0s’*“W 
(1  +  Cl)  Zg 


(11-80) 


For  a  given  gas  source  rate  E,  c  (x)  is  determined  once  values  of 
Sy ,  Sz,  and  b  are  known. 


III.  FORMULATION  OF  THE  DENSE  GAS  DISPERSION  (DEGADIS) 
MODEL  FOR  THE  U.S.  COAST  GUARD  HAZARD  ASSESSMENT 
COMPUTER  SYSTEM  (HACS) 


The  DEGADIS  (DEnse  GAs  Dispersion)  model  is  an  adaptation  of 
the  Shell  HEGADAS  model  described  by  Colenbrander  (1980)  and 
Colenbrander  and  Puttock  (1983).  DEGADIS  also  incorporates  some 
techniques  used  by  van  Ulden  (1983).  If  the  primary  source 
release  rate  exceeds  the  maximum  atmospheric  takeup  rate,  a  HTAG 
blanket  is  formed  over  the  primary  source.  This  near-field, 
buoyancy-dominated  regime  is  modeled  using  a  lumped  parameter 
model  of  a  HTAG  "secondary  source"  cloud  which  incorporates  air 
entrainment  at  the  gravity-spreading  front  using  a  frontal  entrain¬ 
ment  velocity.  If  the  primary  source  release  rate  does  not  exceed 
the  maximum  atmospheric  takeup  rate,  the  released  gas  is  taken  up 
directly  by  the  atmosphere  and  dispersed  downwind.  For  either 
source  condition,  the  downwind  dispersion  phase  of  the  calculation 
assumes  a  power  law  concentration  distribution  in  the  vertical 
direction  and  a  modified  Gaussian  profile  in  the  horizontal  direction 
with  a  power  law  specification  for  the  wind  profile  (Figure  1 1 1 . 1 ) . 

The  source  model  represents  a  spatially  averaged  concentration  of 
gas  present  over  the  primary  source,  while  the  downwind  dispersion 
phase  of  the  calculation  models  an  ensemble  average  of  the  concen¬ 
tration  downwind  of  the  source. 

III.l  HTAG  Source  Cloud  Formation 

A  lumped  parameter  model  of  the  formation  of  the  HTAG  source 
cloud  or  blanket,  which  may  be  formed  from  a  primary  source  such  as 
an  evaporating  liquid  pool  or  otherwise  specified  ground  level 
emission  source,  or  by  an  initially  specified  gas  volume  of  prescribed 
dimensions  for  an  instantaneous  release,  is  illustrated  in  Figure 
III.l.  The  gas  blanket  is  represented  as  a  cylindrical  gas  volume 


Figure  III.l.  Schematic  diagram  of  HACS-DEGADIS 
HTAG  dispersion  model. 
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which  spreads  laterally  as  a  density-driven  flow  with  entrainment 
from  the  top  of  the  source  blanket  by  wind  shear  and  air  entrain¬ 
ment  into  the  advancing  front  edge.  The  source  blanket  will 
continue  to  grow  over  the  primary  source  until  the  atmospheric 
takeup  rate  from  the  top  is  matched  by  the  air  entrainment  rate 
from  the  side  and,  if  applicable,  by  the  rate  of  gas  addition  from 
under  the  blanket.  Of  course,  the  blanket  is  not  formed  if  the 
atmospheric  takeup  rate  is  greater  than  the  evolution  rate  of  the 
primary  source.  For  application  of  the  downwind  calculation 
procedure,  the  blanket  is  modeled  as  being  stationary  over  the 
center  of  the  source  (x  =  0). 

I I I. 1.1  Extent 

If  a  HTAG  blanket  is  present,  the  emission  rate  from  the 
blanket  is  the  maximum  atmospheric  takeup  rate.  That  is,  for 
E(t)  /iTRp(t)  >  Q*,^*  a  source  blanket  is  formed  over  the  liquid 
source.  The  blanket  frontal  (spreading)  velocity  is  modeled  as 


where  p  is  the  average  density  of  the  source  blanket  and  is 
applicable  only  for  p  >  p_;  the  value  of  Cc  used  is  1.15,  based 
on  measurements  of  cloud  spreading  velocity  (Volume  II).  This 
gravity  intrusion  relationship  will  overpredict  velocities  for 
instantaneous,  aboveground  releases  since  no  initial  acceleration 
phase  is  included.  In  this  case,  the  following  procedure  adapted 
from  van  Ulden  (1983)  is  recommended. 

For  instantaneous  gas  releases,  the  radially  symmetric  cloud 
is  considered  to  be  composed  of  a  tail  section  with  height  Ht  and 
radius  Rh  and  a  head  section  with  height  Hh»  A  momentum  balance 
is  used  to  account  for  the  acceleration  of  the  cloud  from  rest 


exhibited  for  experiments  such  as  this  work  (Volume  II)  and  the 
British  Health  and  Safety  Executive's  Thorney  Island  Trials 
(Section  IV. 6)  described  by  McQuaid  (1983).  Since  this  area  of  the 
model  is  meant  to  describe  the  near-field  buoyancy-dominated  flow 
regime,  the  effect  of  any  ambient  momentum  will  be  ignored.  There 
are  three  main  forces  acting  on  the  cloud  including:  a  static 
pressure  force  (Fp),  a  dynamic  drag  force  (F^),  and  a  force  which 
accounts  for  the  acceleration  reaction  of  the  ambient  fluid 
represented  as  a  rate  of  virtual  momentum  change  with  respect  to 
r  dP 

time  -  .  Denoting  the  momentum  of  the  head  and  tail  as 

p,  and  P.  respectively,  the  momentum  balance  becomes 
h  t 


The  terms  in  the  momentum  balance  are  evaluated  differently  for 
early  times  before  a  gravity  current  head  has  developed  and  for 
times  after  the  head  has  developed  (Figure  II. 2)  as  described  in 
Section  II. 1.2,  Equations  11-44  -  65. 

For  early  times,  some  earlier  simplifying  assumptions  do  not 
apply.  In  order  to  model  the  initial  cloud  shape,  the  tail  and 
head  height  are  considered  constant  with  respect  to  radius.  The 
momentum  balance  on  the  cloud  is  given  by 


<Jt  [ph  +  pt]  = 


*9<p  -  “a1  CRh  Ht  +  W  HhJ 


"  TravdvDaR4h  uf 


(III-2) 


where  the  first  term  on  the  right-hand  side  represents  the  static 
pressure  force  on  the  head  and  the  second  term  represents  the  drag 
force  on  the  bottom  surface  of  the  cloud.  A  third  force  is  an 
acceleration  reaction  by  the  ambient  fluid  which  is  represented  by 


dPv  d(RH2uf) 

~dt"  =  ev*pa  dt 
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Figure  III. 2.  Schematic  diagram  of  a  radially  spreading 
cloud. 

The  dimensions  of  the  head  are  determined  by 


Rh*R- 


avbv  Hh 


( 1 1 1 -4 ) 


and 


[g(p  -  p.)  / p - 3 


(III- 


The  constants  ay,  by,  dy,  and  ey  are  assigned  values  1.3,  1.2,  0.64, 
and  20.,  respectively,  based  on  analysis  of  the  still-air  HTA6 
release  experiments  described  in  Section  VI.  When  the  height  of  the 
tail  H+  is  assumed  uniform  with- respect  to  radius,  it  follows  that 

t  a 


Ht'  "avbv(R+  Rh)»h  /(’Rh> 


(III -6 ) 
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where  M  is  the  total  mass  of  the  cloud.  The  momentum  of  the  head 
p^  and  tail  P^  are  determined  by 


.  .2  .  °VR3-Rh> 
’h-rav - : - 


(III -7) 


2  o  Hf 

pt =  r— uf 


(III -8 ) 


Equations  I I I -2  through  II 1-8  determine  the  momentum  of  the 
blanket  as  a  function  of  time,  and  thus  the  frontal  velocity  uf. 
The  cloud  accelerates  from  rest  because  =  0  initially;  when  the 
cloud  accelerates  to  the  point  that  Hh  >  Ht,  the  frontal  velocity 
is  determined  from  Equations  11-44  through  11-65  using  the  values 
of  the  constants  ay,  by,  dy,  and  ey  given  above.  When  the 
frontal  velocity  from  the  momentum  balance  is  the  same  as  Equation 
1 1 1 -1 ,  the  momentum  balance  is  no  longer  applied  and  the  frontal 
velocity  is  given  by  Equation  1 1 1 - 1 . 

The  blanket  radius  R  as  a  function  of  time  is  determined  by 


(III-9) 


dR  =  u 
dt  uf 


For  ground  level  sources,  the  blanket  spreading  is  stopped  ^  *  0 
when  the  total  mass  of  the  cloud  is  decreasing  with  time.  The 
radius  of  the  blanket  is  restricted  to  be  greater  than  or  equal  to 
the  radius  of  the  primary  (liquid)  source  Rp  if  it  is  present. 


I I I. 1.2  Material  and  Energy  Balances 


The  balance  on  the  total  mass  of  gas  in  the  source  blanket 


(M  =  irR2Hp)  is 


$  ■  PHp]  *  «t)  +  +  Aw,s  -  “V2-  H2]  (III'10) 


where  E(t)  is  the  gas  evolution  rate  from  the  primary  (liquid) 
source.  For  spills  over  water,  the  water  entrainment  term  (M  ) 

W  »  5 

is  included  in  the  source  blanket  description  and  is  calculated  from 
Equation  III -24 ,  and  the  air  entrainment  rate  is 


Mfl  ■  2  tt  RHh(euf)  pa  j  [g  A'  H  /  u2] 


(III-ll) 


the  value  of  e  is  estimated  to  be  0.6,  based  on  the  still-air  HTAG 
release  experimental  data  of  Volume  II. 

The  balance  on  the  mass  of  contaminant  in  the  source  blanket 
(Mp  *  wc  if  r2h p)  is 


dM 

ir-sfc’"2"']  -E<‘)  -<W*r2> 


( I I I -1 2) 


and  the  mass  balance  on  the  air  in  the  source  blanket  (M^  *  w.itITHp) 
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where  the  ambient  humidity  is  H.  and  the  mass  fraction  of  contaminant 

d 

and  air  are  w„  s  M„  /  M  and  =  Ma/M,  respectively. 

C  C  a  a  n 

The  energy  balance  on  the  source  blanket  (hirRcHp)  gives 


'SSmsmsmmaiBiSm 


60 


^  h7rR2HpJ  = 


hEE(t) 


haMa 


hwMw,s 


-  h 
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(III-14) 


where  h£  is  the  enthalpy  of  the  emitted  gas,  ha  is  the  enthalpy  of 
the  ambient  humid  air,  and  hw  is  the  enthalpy  of  any  water  vapor 
entrained  by  the  blanket  if  over  water.  There  are  three  alternate 
submodels  included  for  the  heat  transfer  (Qs)  from  the  surface  to 
the  cloud. 

The  simplest  method  for  calculating  the  heat  transfer  between 
the  substrate  and  the  gas  cloud  is  to  specify  a  constant  heat 
transfer  coefficient  for  the  heat  transfer  relation 


Qs  *  qs  f(R2  -  R2)]  *  h04T  |fl(R2  -  R2)]  (III-15) 

where  Qs  is  the  rate  of  heat  transfer  to  the  cloud,  q$  is  the  heat 
flux,  and  AT  is  the  temperature  difference.  For  the  calculation  of 
heat  transfer  over  the  source,  the  temperature  difference  is  based 
on  the  average  temperature  of  the  blanket. 

In  the  evaluation  of  the  Burro  and  Coyote  series  of  experiments, 
Koopman  et  al.  (1981)  proposed  the  following  empirical  heat  transfer 
coefficient  relationship  for  heat  transfer  between  a  cold  LNG  cloud 
and  the  ground 

*0*VCp  (III-16) 

where  the  value  of  was  estimated  to  be  0.0125  m/s.  This  constant 
can  be  varied  in  the  model. 

From  the  heat  transfer  coefficient  descriptions  for  heat  transfer 
from  a  flat  plate,  the  following  relationships  can  be  applied.  For 


natural  convection,  the  heat  transfer  coefficient  is  estimated 
using  the  Nusselt  Nu,  Grashoff  Gr,  and  Schmidt  Sc  numbers,  McAdams 
(1954),  from 


0.14  (GrSc)1/3 


(111-17) 


r  2  3 

9  p  C  y 

h  =  0.14  - X-  AT 

n  TPr2 


( II I -18) 


where  h  is  the  heat  transfer  coefficient  due  to  natural  convection 
n 


and  Pr  is  the  Prandtl  number.  In  order  to  simplify  the  calculations, 
the  parameter  group 


Or*2  (Cp  MW)3  (pjfr) 


(III— 19) 


is  estimated  to  be  60  in  mks  units.  The  actual  value  of  the  group 
is  47.25,  58.5,  and  73.4  for  air,  methane,  and  propane,  respectively. 
Equation  I 11-18  then  becomes 


h„*  '8 


(III -20) 


where  the  density  p,  molecular  weight  MW,  and  temperature  difference 
AT  are  based  on  the  average  composition  of  the  gas  blanket. 

For  forced  convection  after  Treybal  (1980),  the  Colburn  analogy 
is  applied  to  a  flat  plate  using  the  Stanton  number  for  heat  transfer 
StH  and  the  Prandtl  number  as 


St.  Pr2^3  *  —  -  h\ 


(III-21 ) 
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Or, 


if  a  ( JoCp)  Pr 
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(III -22 ) 


where  h^  is  the  heat  transfer  coefficient  due  to  forced  convection. 

If  the  velocity  is  evaluated  at  the  top  of  the  gas  blanket  and  the  Pr 
is  estimated  to  be  0.741, 
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The  overall  heat  transfer  coefficient  is  then  the  maximum  of  the 
forced  and  natural  coefficients,  i.e.  hQ  =  max(hf,  hn).  The  heat  flux 
and  transfer  raie  are  then  estimated  by  Equation  III-15. 

If  the  gas  blanket  is  formed  over  water,  water  will  be  transferred 
from  the  surface  to  the  cloud  by  a  partial  pressure  driving  force 
associated  with  the  temperature  difference  between  the  surface  and 
the  gas  blanket.  The  rate  of  mass  transfer  of  water  is 


’O2  -  Rp)J 


(III -24 ) 


where  Fq  is  the  overall  mass  transfer  coefficient.  The  driving  force 

is  the  difference  of  the  vapor  pressure  of  water  at  the  surface 

temperature  p*  and  the  vapor  pressure  of  water  at  the  cloud  tempera- 
*  w»s 

ture  p  .  The  natural  convection  coefficient  is  based  on  the  heat 
w,c 

transfer  coefficient  and  the  analogy  between  the  Sherwood  number  Sh 
and  the  Nusselt  number  Nu  suggested  by  Bird  et  al .  (1960);  or 


Sh  *  Nu  *  0.14  (GrSc)1/3 


(III -25) 


where  Sc  is  the  Schmidt  number.  If  the  Schmidt  number  is  taken  as 
0.6,  and 
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T  MW 


is  estimated  to  be  2.2  x  10"®  in  mks  units. 


Fn  »  9.9  x  10 


-3 


itf  *t 


-rl/3 


(III -26 ) 


For  forced  convection,  Treybal  (1980)  suggests  that  the  Stanton  number 
for  mass  tr 
related  by 


for  mass  transfer  StM  and  the  Stanton  number  for  heat  transfer  StH  are 


StM  *  StH 
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Or, 
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(III —28 ) 


The  overall  ma  s  transfer  coefficient  Fq  is  calculated  as  the  larger 
of  the  natural  and  forced  convection  coefficients. 

For  the  c< se  when  the  primary  (liquid)  source  emission  rate  E(t) 

.  2 

is  larger  than  the  atmospheric  takeup  rate  Q+max^Rp*  Equations 

III-10,  1 1 1 -1 2  III-13,  and  1 1 1 -1 4  are  integrated  for  the  mass, 
concentration,  and  enthalpy  of  the  gas  blanket  along  with  an 
appropriate  eqi  ation  of  state  (i.e.  relationship  between  enthalpy  and 
temperature  an>  between  temperature  and  density)  from  Appendix  A. 

For  the  cjse  when  the  emission  rate  is  insufficient  to  form  a  gas 
blanket,  the  flux  of  contaminant  is  no  longer  determined  by  the 
maximum  atmospheric  takeup  rate.  Consider  the  boundary  layer  formed 
by  the  emission  of  gas  into  the  atmosphere  above  the  primary  source. 

If  the  source  is  modeled  to  have  a  uniform  width  2b  and  entrain  no 
air  along  the  sides  of  the  layer,  the  balance  on  the  total  material 
(pL  uLHL^  a  d^'ferent'ia1  slice  of  the  layer  is 
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(III -29) 


where  wg  is  the  vertical  rate  of  air  entrainment 
by  Equation  III -61 »  is  the  average  density  of 
Q*  /  wc  is  the  total  flux  of  gas  from  the  primary 
The  balance  on  the  mass  flow  rate  of  contaminant 


into  the  layer  given 
the  slice,  and 
(liquid)  source. 
(wcplulHl)  at  any 
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With  an  equation  of  state  to  relate  c  ,  and  p.  ,  Equation  1 1 1 -29 

u  )  L  L  . 

is  integrated  from  the  upwind  edge  of  the  source  (x  =  xup)  to 
the  downwind  edge  (x  =  L  +  xup). 

In  order  to  generate  the  initial  conditions  for  the  downwind 
dispersion  calculations,  the  maximum  concentration  cc  and  the  vertical 
dispersion  parameter  Sz  are  needed.  Since  Equations  1 1 1-29  and  II 1-30 
are  written  for  a  vertically  averaged  layer,  consider  the  vertical 
average  of  the  power  law  distribution.  The  height  of  the  layer 
is  the  height  to  some  concentration  level,  say  10%  of  the  maximum. 
Although  strictly  a  function  of  a,  this  value  is  modeled  by 


Hl  *  6lHeff  ( 1 1 1 -31 ) 

where  HEFF  is  the  effective  height  defined  by  Equation  1 1 1-57  and 
6l  is  2.15.  The  vertically  averaged  concentration  cc  L  can  be 
defined  by 
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c,L 


(III -32 ) 


And  similarly,  the  effective  transport  velocity  is  defined  by 


OD 

:c,iuiA  ’f  cu> 


With  Equation  1 1 1-31  and  defining  relations  for  and 
(Equations  III -57  and  II I -71,  respectively),  it  follows  that 
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where  w'  is  given  by  Equation  I I 1-61 


I I I. 1.3  Maximum  Atmospheric  Takeup  Rate 


(III-- 


(in-: 


The  maximum  atmospheric  takeup  rate  will  be  the  largest  takeup 
rate  which  satisfies  Equations  III -29  and  III-30.  As  well,  the 
maximum  concentration  of  contaminant  in  the  power  law  profile  at  the 
downwind  edge  of  the  source  will  be  the  source  contaminant  concentra 
tion  (cc)s-  If  Equations  III -29  and  1 1 1-30  are  combined  along  with 
the  assumption  of  adiabatic  mixing  of  ideal  gases  with  the  same 


molal  heat  capacity  (i.e.  — 


p  -  p. 


*  Ys  constant),  the  maximum 


takeup  flux  is  modeled  by 


k  u*  (1  +  a)  f  <5,  1 
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where 


L 

1  _  1  f  dx 

♦  '  L  J  * 


(III-38) 


An  upper  bound  of  the  atmospheric  takeup  flux  can  be  characterized 
by  the  condition  where  the  source  begins  to  spread  as  a  gravity 
intrusion  against  the  approach  flow.  In  water  flume  experiments, 
Britter  (1980)  measured  the  upstream  and  lateral  extent  of  a  steady 
state  plume  from  a  circular  source  as  a  function  of  Ri*.  A 
significant  upstream  spread  was  obtained  for  Ri*  >  32,  and  lateral 
spreading  at  the  center  of  the  source  was  insignificant  for  Ri*  <  8. 

The  presence  of  any  significant  lateral  spreading  represents  a  lower 
bound  on  the  conditions  of  the  maximum  takeup  flux. 

The  integral  of  Equation  III -38  is  calculated  using  a  local 
Richardson  number  of 


Ri*  (x)  *  c(x  -  xup) 


(III -39) 


where 


’TTJ  ku*<1  +  a>fw  ](  SL 


1  +  a 


u0  Z0 


(III -40) 


and  where  the  value  of  <pc  is  3.1  (corresponding  to  Ri*  =  20 
(8  <  Ri*  <32)).  Using  this  Ri  *  ( x )  and  the  first  two  terms  of 
4>(Ri*),  Equation  1 1 1 -38  is  approximated  as 


<D  0.099  L  £ 


In  0-88  + 0.099c 

0.88 


(III -41 ) 


mm 


which  then  specifies  the  maximum  atmospheric  takeup  flux. 


III. 1.4  Transient  HTAG  Release  Simulation 


If  a  steady  state  spill  is  being  simulated,  the  transient  source 
calculation  is  carried  out  until  the  source  characteristics  are  no 
longer  varying  significantly  with  time.  The  maximum  centerline 
concentration  cc,  the  horizontal  and  vertical  dispersion  parameters 
$y  and  S2,  the  half  width  b,  and  if  necessary,  the  enthalpy  h  are 
used  as  initial  conditions  for  the  downwind  calculation  specified  in 
a  transient  spill. 

If  a  transient  spill  is  being  simulated,  the  spill  is  modeled 
as  a  series  of  pseudo-steady  state  releases.  Consider  a  series  of 
observers  traveling  with  the  wind  over  the  transient  gas  source 
described  above;  each  observer  originates  from  the  point  which 
corresponds  with  the  maximum  upwind  extent  of  the  gas  blanket 
(x  =  ’Rmax)*  The  desired  observer  velocity  is  the  average  transport 
velocity  of  the  gas  u^p  from  Equation  III -71 ;  however,  the  value  of 
UEFF  w1^  differ  from  observer  to  observer  with  the  consequence  that 
some  observers  may  be  overtaken  by  others.  For  a  neutrally  buoyant 
cloud,  Uppp  becomes  a  function  of  downwind  distance  alone  which 
circumvents  this  problem.  With  this  functionality,  Colenbrander 
(1980)  models  the  observer  velocity  as 


u^x) 


r+- 


'0m 


a 


x  +  R. 


-ret  /  (1+a) 


max 


^"Rm  +  Rmax 


(III -42 ) 


•  ,  ? 
where  S,  is  the  value  of  S  when  the  averaged  source  rate  '.it  R  Q*) 

z0m  z0 

is  a  maximum, and  the  subscript  i  denotes  observer  i.  Noting  that 
dx^. 

Ui(x)  =  -g^r  .  observer  position  and  velocity  as  functions  of  time  are 
determined. 
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A  pseudo-steady  state  approximation  of  the  transient  source  is 

obtained  as  each  observer  passes  over  the  source.  If  t  and  t 

up,.  dn^ 

denote  the  times  when  observer  i  encounters  the  upwind  and  downwind 
edges  of  the  source  respectively,  then  the  source  fetch  seen  by 
observer  i  is: 


Li  *  xup.  "  xdn. 


(Ill -43) 


The  width  of  the  source  2Bj(t)  is  defined  by 


sf(t)  =■  R2(t)  -  x2(t) 


(III -44 ) 


Then,  the  gas  source  area  seen  by  observer  i  is 

t. 


dn,. 


2Libi  *  2  /  Bj  «,  dt 


(III -45 ) 


'up,- 


where  2b..  is  the  average  width. 

The  takeup  rate  of  contaminant  (Q*Lb)i  is  calculated  by 

t. 


dn. 


2(Q 


.Lb),  ■  2  J 


Q*  Bi  u.  dt 


(I 11-46) 


'UP; 


The  total  mass  flux  rate  from  the  source  is 


dni 


2(pLuLHLb).  =2 


/ 


pawe 


Q* 


Bi  dt 


(III -47 ) 
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With  these  equations,  the  average  composition  of  the  layer  can 
be  determined  at  each  x  *  xup  over  the  source.  With  the  enthalpy 
of  the  layer  given  by 


(hLPLuLHLb), 


(in-48) 


due  to  the  choice  of  the  reference  temperature  as  the  ambient 
temperature  along  with  a  suitable  equation  of  state  relating 
enthalpy,  temperature,  and  density,  the  source  can  be  averaged  for 
each  observer.  After  the  average  composition  of  the  layer  is 
determined  at  the  downwind  edge,  an  adiabatic  mixing  calculation 
is  performed  between  this  gas  and  the  ambient  air  when  applicable. 

This  calculation  represents  the  function  between  density  and 
concentration  for  the  remainder  of  the  calculation  if  the  calculation 
is  adiabatic;  it  represents  the  adiabatic  mixing  condition  if  heat 
transfer  is  Included  in  the  downwind  calculation. 

For  each  of  several  observers  released  successively  from 
x  3  -R_v,  the  observed  dimensions  L  and  b,  the  downwind  edge  of 
the  source  xdn,  the  average  vertical  dispersion  coefficient  Sz,  the 
average  takeup  flux  Q*,  the  centerline  concentration  cc,  and  if 
applicable,  the  average  enthalpy  hL  can  be  determined  for  each 
observer.  With  these  input  values,  a  steady  state  calculation  is 
made  for  each  observer.  The  distribution  parameters  for  any  specified 
time  t$  are  determined  by  locating  the  position  of  the  series  of 
observers  at  time  t$,  i.e.  x^(t$).  The  corresponding  concentration 
distribution  is  then  computed  from  the  assumed  profiles. 
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I I I. 2  Steady  State  Downwind  Dispersion 

The  model  treats  dispersion  of  gas  entrained  into  the  wind  field 
from  an  idealized  rectangularly  shaped  source  of  width  2b  and  length 
L.  The  circular  source  cloud  is  represented  as  an  equivalent  area 
square  (L  *  ttR  *  2bL).  Similarity  forms  for  the  concentration 
profiles  are  assumed  which  represent  the  plume  as  being  composed  of 
a  horizontally  homogeneous  section  in  which  only  vertical  dispersion 
is  considered,  with  Gaussian  concentration  profile  edges  as  follows: 


c(x,y,z)  =  cc(x)exp 


|y|  -  b(x) 

2 

'  \ 

z 

1+a 

Sy(x) 

k  ~  • 

Sz(x) 

for  jy |  >  b 


cc(x)  exp 


*  ' 

2 

l+a 

»  * 

for  |y  |  js  b 
(III -49) 


A  power  law  wind  velocity  profile  is  assumed 


U  a  Un 
x  0 


(III -50) 


where  the  value  of  a  is  determined  by  a  weighted  least  squares  fit 
of  the  logarithmic  profile 


In 


z  +  zf 


:R 

I 


* 


r  \ 

i 

X 
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Functional  forms  for  ip  and  typical  values  of  a  are  given  in  Table 
II I -1  for  different  Pasquill  stability  categories.  With  these 
profiles,  the  parameters  of  Equation  II 1-49  are  constrained  by 
ordinary  differential  equations. 

I I I. 2.1  Vertical  Dispersion 


The  vertical  dispersion  parameter  S2  is  determined  by  requiring 
that  it  satisfy  the  diffusion  equation 


3c 


ii  IS.  *  ^  k 

ux  3x  37  KzSz 


(III -52 ) 


with  the  vertical  turbulent  diffusivity  given  by 


ku*z 


amj 


(III -53) 


The  function  $(Ri*)  is  a  curve  fit  of  laboratory  scale  data 
for  vertical  mixing  in  stably  density-stratified  fluids  reported  by 
McQuai^  (1976),  Kantha  et  al.  (1977),  and  Lofquist  (1960)  for 
R1*  >  0.  For  R1*  <0,  the  function  <|>(R1*)  is  taken  from  Colenbrander 
(1983)  and  has  been  modified  so  the  passive  limit  of  the  two  functions 
agree  as  follows: 


<fc(R1*)  *  0.88  +  0.099  Ril’04  +  1.4  x  10’25  Ri*'7  Ri*  >  0 


0.38/  (1  +  0.65  | Ri*|0'6) 


Ri*  <  0 

(III -54 ) 


The  friction  velocity  is  calculated  using  Equation  III -51  from  a 
known  velocity  uQ  at  a  specific  height  zQ.  Combining  the  assumed 
similarity  forms  for  concentration  and  velocity.  Equations  III -49 , 

II 1-50,  1 1 1-52,  and  III-53  give 


Curve  fit  of  data  from  Pasquill  (1974). 
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d  >Q  Zq  fSz]l4<J|  k  u*  (1  +  a) 
dx  r+"£ 

where  the  Richardson  number  Ri*  is  computed  as 


(III -55 ) 


Ri*  *  g 


p  -  P,  H, 


a  EFF 


a  u* 


(II 1-56 ) 


and  the  effective  cloud  depth  is  defined  as 


57  jcdz  ’ r  (t^t)  rr 


(III -57 ) 


1!  nR 

JJS 


Equation  111-55  can  be  viewed  as  a  volumetric  balance  on  a 
differential  slice  of  material  downwind  of  the  source.  For  a  mass 
balance  over  the  same  slice,  one  obtains 


3x  IpLulhl  pa 


(III -58 ) 


which  is  the  same  result  as  Equation  I 11-29  without  the  source  term. 
With  Equations  II 1-35  and  III-36,  this  is 


dx  pLuEFFHEFF  *  pa"e 


(III -59 ) 


Using  the  assumption  of  adiabatic  mixing  of  ideal  gases  with  the 

p  -  p 

same  constant  molal  heat  capacity  (i.e.  — - - =  constant)  along  with 

the  contaminant  material  balance,  the  mass  balance  becomes 
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d 

37 


UEFF  HEFF 


which  leads  to 


(I I 1-60) 


(III -61 ) 


Equations  I I 1-59  and  I 11-61  are  combined  to  give 


d 

37 


Paku*(l  +  a) 


(III -62 ) 


P  -  P. 

for  situations  even  when  — - — -  is  not  constant. 

cc 

When  heat  transfer  from  the  surface  is  present,  vertical  mixing 
will  be  enhanced  by  fc.he  convective  turbulence  due  to  heat  transfer. 
Zeman  and  Tennekes  (1977)  model  the  resultant  vertical  turbulence 
velocity  (Equation  I I -5)  as 
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' 

w* 

3 

k  4 

-|l/2 


where  w*  is  the  convective  scale  velocity  described  by 

—r  2/3 


’** 

q  h  ^s  Tc,lJ 

k  . 

u*u  Tc,L 

If  u  is  evaluated  at  HEFF’ 
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where 


RiT  *  9 


Ts  "  Tc.L  HEFF  20 
T  .  U*U0  HEFF 
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and  T_  ,  is  the  temperature  obtained  from  the  energy  balance  of 
Equations  III -80  and  I II -81.  Equation  I I 1-62  is  modified  to  account 
for  this  enhanced  mixing  by 


a*  |pLuEFFHEFF 


Pa  kw(l  +  a) 

d>(Ri*) 


(II 1-67) 


where  Ri*  =  RiJ— j  * 


Although  derived  for  two-dimensional  dispersion,  this  is  extended 
for  application  to  a  HTA6  plume  which  spreads  laterally  as  a  density 
intrusion: 


2x  °L  UEFF  HEFF  BEFF 


p  k  w(l  +  a) 

-J - B, 
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(III -68) 


where  the  plume  effective  half  width  is  defined  by 


BEFF  *  b  +  T*  Sy 


(III -69) 


and  determined  using  the  gravity  intrusion  relation 


cf  9 


P  -  P, 
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The  average  transport  velocity  in  the  plume  is  defined  by 


OD 
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f  cdz 
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and  the  lateral  spread  of  the  cloud  is  modeled  by 


dBEFF  .  1  d  BEFF 

dx  UEFF  dt 


1/2  fsJ(*-a) 
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(II 1-72 ) 


I I I. 2. 2  Horizontal  Dispersion 


The  crosswind  similarity  parameter  S^(x)  is  also  determined  by 
requiring  that  it  satisfy  the  diffusion  equation 


3c  -  _L  «  i£ 
ux  3x  ”  3y  Ny  3y 


with  the  horizontal  turbulent  diffusivity  given  by 


Ky  3  K0  ux  BEFF 
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When  b  =  0,  Su  ■  /T  o  ,  where  a  is  the  similarity  parameter 

y  *  *  e 

correlated  by  Pasquill  (1974)  in  the  form  =  6x  .  Furthermore, 
Equations  1 1 1-73  and  1 11-74  require  that 


Jy  ■  "o  beff 


(I  II -75) 


where  y-j  =2-1/8  and  Kg  *  (<5/  ir/2  )  Then 


,  th.  *e  B  2  {/i7r 
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where  Equation  III-76  is  also  assumed  applicable  for  determining 
Sy  when  b  is  not  zero. 

At  the  downwind  distance  x^  where  b  s  0,  the  crosswind 
concentration  profile  is  assumed  Gaussian  with  Sy  given  by 


Sy  *  /T  6 ( x  +  xv)f 


(II 1-77 ) 


where  xy  is  a  virtual  source  distance  determined  by 


Sy(xt)  »  /T  6(xt  ♦  x/ 


( 1 11-78) 


The  gravity  spreading  calculation  is  terminated  for  x  >  xt> 

For  a  steady  plume,  the  centerline  concentration  cc  is  determined 
from  the  material  balance 
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where  E  is  the  plume  source  strength. 
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I I I. 2. 3  Energy  Balance 

For  some  simulations  of  cryogenic  gas  releases,  heat  transfer 
to  the  plume  in  the  downwind  dispersion  calculation  may  be  important, 
particularly  in  low  wind  conditions.  The  source  calculation 
determines  a  gas/air  mixture  initial  condition  for  the  downwind 
dispersion  problem.  Air  entrained  into  the  plume  is  assumed  to 
mix  adiabatically.  Heat  transfer  to  the  plume  downwind  of  the  source 
adds  additional  heat.  This  added  heat  per  unit  mass  Dh  is  determined 
by  an  energy  balance  on  a  uniform  cross-section  as 


lx  Dh  PL  UEFF  HEFF  *  qs  /  5l 


(I I 1-80) 


where  qs  is  determined  by  Equation  I 11-15  along  with  the  desired 
method  of  calculating  hQ.  Equation  I I 1-80  is  applied  when  b  *  0  and 
is  extended  to 


[°h  PL  UEFF  HEFF  BEFfJ  *qsBEFF/6l 


(II 1-81 ) 


when  b  >  0.  Since  the  average  density  of  the  layer  pL  cannot  be 
determined  until  the  temperature  (i.e.  Dh)  is  known,  a  trial  and 
error  procedure  is  required. 

Equations  III -55 ,  56,  57,  65-69,  72,  and  76-81  are  combined  with 

an  equation  of  state  relating  cloud  density  to  gas  concentration  and 

temperature  and  are  solved  simultaneously  to  predict  S  ,  S  ,  cr,  and 

z  y  c 

b  as  functions  of  downwind  distance  beginning  at  the  downwind  edge 
of  the  gas  source. 


EsfflSr&Ci-a 
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I I I. 3  Correction  for  Along-Wind  Dispersion 


In  a  similar  fashion  to  Colenbrander  (1980),  an  adjustment  to 
the  values  of  cc  is  applied  to  account  for  dispersion  parallel  to 


the  wind  direction.  The  calculated  centerline  concentration  cc (x) 


is  considered  to  have  resulted  from  the  release  of  successive  planar 
puffs  of  gas  (cc(x)Ax)  without  any  dispersion  in  the  x-direction.  If 
it  is  assumed  that  each  puff  diffuses  in  the  x-direction  as  the  puff 
moves  downwind  independently  of  any  other  puff  and  that  the  dispersion 
is  one-dimensional  and  Gaussian,  the  x-direction  concentration 
dependence  is  given  by 
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where  xp^  denotes  the  position  of  the  puff  center  due  to  observer  i. 

After  Beals  (1971),  the  x-direction  dispersion  coefficient  ax 
is  assumed  to  be  a  function  of  distance  from  the  downwind  edge  of  the 
gas  source  (X  =  x  -  xQ)  and  atmospheric  stability  given  by 


<Jx(X)  *  0.02  X1-22 

unstable. 

■  0 . 04  X 1  ‘ 1 4 
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stable,  . 
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where  (X  =  x  -  Xg)  and  ax  are  in  meters.  The  concentration  at  x  is 
then  determined  by  superposition,  i.e.,  the  contribution  to  c  at  a 
given  x  from  neighboring  puffs  is  added  to  give  an  x-direction 


corrected  value  of  c  .  For  N  observers, 
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and  for  large  N, 
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The  corrected  centerline  concentration  c^  is  used  in  the  assumed 

profiles  in  place  of  cr,  along  with  the  distribution  parameters 

c  y 

S  ,  and  b. 
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III. 4  The  DEGADIS  Computer  Model 

The  model  described  in  the  previous  section  has  been  programmed 
in  Digital  Equipment  Company  (DEC)  VAX  FORTRAN  and  resides  on  a  DEC 
VAX  11-730  under  VMS. 

Volume  III  of  this  report  contains  a  flow  chart  of  the  DEGADIS 
computer  model,  a  complete  FORTRAN  code  listing,  and  an  example 
DEGADIS  interactive  computer  session. 


IV.  COMPARISON  OF  SIMULATIONS  WITH 
SELECTED  FIELD  EXPERIMENT  DATA 


The  DEGADIS  model  has  been  used  to  simulate  a  collection  of 
thirty-nine  field  experimental  HTA6  releases.  The  tests  simulated 
include  small,  continuous  LPG  releases  (order  0.1  -  1.0  kg/s)  on 
land  from  diked  sources,  continuous  LPG  and  LNG  releases  onto  water 
with  release  rates  of  the  order  10  -  100  kg/s,  near-instantaneous 
releases  of  LNG  of  approximately  5000  kg  onto  water,  and 
instantaneous  releases  of  Freon-12/air  mixtures  of  approximately 
5000  kg  on  land.  The  meteorological  conditions  encompassed  wind 
speeds  of  approximately  1-10  m/s,  relative  humidity  from 
essentially  zero  to  about  85%,  and  estimated  atmospheric 
stabilities  ranging  from  Pasquill  B  to  F. 

In  the  following  sections,  the  individual  HTAG  dispersion  test 
series  simulated  are  briefly  described,  and  the  simulation  results 
are  compared  with  field  experimental  measurements.  For  those  tests 
which  can  be  reasonably  represented  as  either  instantaneous  or 
steady  continuous  releases,  a  "release  Richardson  number,"  Rig, 
is  presented.  This  Richardson  number  may  be  useful  for  demarcating 
HTAG  releases  which  exhibit  significant  departures  from  passive 
dispersion  behavior  (Puttock,  1982).  The  final  section  summarizes 
observations  of  the  agreement  between  reported  field  experimental 
gas  concentration  measurements  and  DEGADIS  model  predictions. 

* 

IV. 1  LPG  Releases  on  Land  from  Diked  Area  Sources 

Welker  (1982)  reported  vapor  dispersion  measurements  for 
propane  releases  into  concrete-lined  pits  nominally  1.5,  3,  and 
6  meters  square  in  a  variety  of  wind  and  atmospheric  stability 
conditions.  Measurements  of  gas  concentration  at  ground  level  at 
positions  estimated  to  be  in  the  centerline  of  the  plume  path  were 
reported. 


Table  IV.l  gives  the  estimated  source  strength  and  meteorological 
data  for  ten  of  these  tests  which  have  been  simulated  with  DEGADIS. 

The  reported  measured  gas  concentrations  (taken  at  a  height  above 
ground  of  0.15  m)  were  obtained  with  a  catalytic  type  hydrocarbon 
sensor  having  a  response  time  estimated  to  be  between  1  and  5  s. 
Figures  IV.l  through  IV. 10  show  the  DEGADIS  model -predicted  ground 
level  centerline  gas  concentration  compared  with  the  reported 
maximum  and  time-averaged  gas  concentrations  for  Welker's  tests  275-1, 
281-2,  296-1,  297-1,  297-3,  298-1,  298-2,  302-1,  334-1,  and  334-2 
respectively. 


TABLE  IV-1 


SUMMARY  OF  SOURCE  AND  METEOROLOGICAL  DATA 
FOR  SELECTED  PROPANE  RELEASES  FROM  WELKER 


Run  No. 

Width  (m) 

Propane  Evap.  Rate 
(kg/s) 

Wind  Speed 
iVs’ 

Stability2 

Air  - 

Temperature J 

00 

Relative 

Humidity 

% 

275-1 

1.52 

0.0189 

4.52 

B 

301.5 

28 

281-2 

3.04 

0.057 

3.09 

B 

302.6 

21 

296-1 

3.04 

0.071 

4.29 

D 

294 

40 

297-1 

3.04 

0.058 

4.34 

C 

296 

23 

297-3 

3.04 

0.051 

3.67 

0 

298 

38 

298-1 

3.04 

0.055 

2.01 

C 

294 

23 

298-2 

3.04 

0.044 

2.15 

B 

296 

30 

302-1 

3.04 

0.073 

7.47 

0 

295 

69 

334-1 

6.10 

0.249 

2.37 

0 

-279 

-10 

334-2 

6.10 

0.551 

2.77 

c 

-279 

-10 

*at  -8.0  ft.  (2.44  m)  height 

^asqulll  Category,  estimated  by  Welker 

Surface  temperature  estimated  T  ■  T  ,  ♦  1  -  3°K 
Surface  roughness  estimated  •  0:01  m 


.t  • 


DEGADIS-predlcted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Welker  275-1 . 
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Figure  IV. 2.  DEGADIS-predlcted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Welker  281-2. 
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DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Welker  296-1. 
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Figure  IV. 4.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Welker  297-1. 
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Figure  IV. 5.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Welker  297-3. 
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Figure  IV. 8.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Welker  302-1. 
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Figure  IV. 9.  OEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
J  Welker  334-1. 
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IV. 2  LNG  Releases  on  Land  from  Diked  Area  Sources 


The  American  Gas  Association  sponsored  a  series  of  more  than 
30  LNG  releases  (AGA,  1974)  into  diked  land  areas  up  to  24.4  m 
diameter  in  1973  near  San  Clemente,  California.  Releases  were 
made  into  1.8  m,  6.1  m,  and  24.4  m  diameter  land  areas  surrounded 
by  ~0.5  m  high,  insulated  wall  dikes.  Vapor  release  rates  for  the 
experiments  were  determined  from  correlations  of  experimentally 
measured  liquid  boiling  rates  in  the  test  program,  and  the  LNG 
vapor  source  rate  used  for  1.8,  6.1,  and  24.4  m  diked  areas  are 
shown  in  Figure  IV. 11.  The  measured  maximum  gas  concentration  (at 
a  height  of  0.5  m)  from  eleven  of  the  experiments,  which  included 
releases  in  all  three  dike  sizes,  were  correlated  with  wind 
velocity,  diked  area,  and  downwind  distance  as  shown  in  Figure 
IV. 12  (AGA,  1974). 

N. 

The  DEGADIS  model  has  been  used  to  simulate  LNG  releases' with 
the  source  rates  and  dimensions  shown  in  Figure  IV. 11  and  with  a 
wind  velocity  of  3.6  m/s,  representing  "typical"  releases  from 
1.8  m,  6.1  m,  and  24.4  m  dikes  respectively.  Atmospheric  stability 
C  was  assumed  with  ambient  and  soil  surface  temperatures  288K  and 
300K,  respectively.  The  ambient  humidity  was  assumed  45%.  The 
predicted  downwind  distance  to  the  LFL  concentration  (5%  at  ground 
level)  for  the  three  simulations  is  shown  in  Figure  IV. 12  for 
comparison  with  the  experimental  data  correlation. 

Figure  IV. 13  shows  the  DEGADIS-predicted  maximum  concentration 
decay  (at  ground  level)  with  distance  from  release  center  compared 
with  reported  experimental  measurements  from  AGA  Test  44,  the 
largest  LNG  release  of  the  test  series,  for  which  the  pertinent 
test  description  is  given  in  Table  IV. 2. 


■.■f.V.V.V. 


DEGADIS  predictions  for  1.8,  6.1,  and  24.4  meter  diked  land  releases  compared 
with  Battel le  correlation  of  maximum  measured  concentrations. 


TABLE  IV. 2 

TEST  CONDITIONS  FOR  AGA  TEST  44 

LNG  Volume  Spilled:  51  m3  (liq)  in  20-30  s 
Diked  Area  Diameter:  24.4  m 

Vapor  Release  Rate:  est.  as  shown  in  Figure  IV. 11 
Wind  Speed:  3.6  m/s  (actual  varied  3.6  -  7.2  m/s) 
Ambient  Temperature:  14.5  C 
Surface  Temperature:  estimated  14.5  C 
Atmospheric  Stability:  estimated  Pasquill  C 
Surface  Roughness:  estimated  1  cm 
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DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
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IV. 3  American  Petroleum  Institute/ESSO  LNG  Releases  on  Water 

The  American  Petroleum  Institute  sponsored  a  series  of  17  LN6 

releases  conducted  by  ESSO  Research  and  Engineering  Company  on 

Matagorda  Bay,  Texas,  in  1972  (API,  1972).  The  releases  included 

3  3 

liquid  volumes  ranging  from  about  1  m  to  10.2  m  with  release 

3  3 

rates  ranging  from  about  10  m  /min  to  about  20  m  /min.  The  LNG 

3 

was  discharged  from  a  nitrogen-pressurized  ~13  m  tank  as  a  jet. 

The  liquid  exit  point  was  elevated  and  about  30°  above  horizontal. 

It  has  been  pointed  out  by  Shell  Research  that  a  substantial  amount 

of  the  LNG  could  have  vaporized  before  hitting  the  sea  surface; 

estimates  of  the  amount  vaporized  in  the  jet  range  as  high  as  75%. 

Since  such  vaporization  would  result  in  an  LNG  vapor/air  mixture 

denser  than  pure  LNG  vapor,  the  description  of  the  source  as  input 

to  a  vapor  dispersion  model  is  uncertain.  Predictions  for  this 

spill  using  DEGADIS  are  consistent  with  the  experimental  results 

when  the  source  is  represented  as  an  instantaneous  spill  of  liquid 

2 

on  water  with  an  LNG  boiling  rate  of  0.196  kg/m  s.  This  boiling 
rate,  which  is  about  twice  as  high  as  indicated  by  measurements  by 
Shell  (Blackmore,  1982),  may  be  consistent  with  the  experiments 
because  it  more  closely  represents  the  rapid  cloud  formation  by  the 
jet. 

Table  IV. 3  shows  the  pertinent  test  descriptions  for  tests  11 

and  17,  which  were  the  two  largest  LNG  releases,  simulated  with  the 

DEGADIS  model.  ESSO  11  and  ESSO  17  were  simulated  as  instantaneous 

liquid  releases  onto  water,  and  the  LNG  liquid  pool  growth  was 

estimated  using  the  liquid  spread  model  described  by  Havens  (1982), 

2 

which  incorporated  an  LNG  boiling  rate  of  0.196  kg/m  s. 

Figures  IV. 14  and  IV. 15  show  the  DEGADIS  model -predicted 
maximum  ground  level  centerline  gas  concentration  compared  with 
measured  maximum  concentrations  (at  ~0.4  m  height)  as  a  function 
of  downwind  distance  from  the  release  center  for  ESSO  11  and  ESSO  17 
respectively.  The  gas  concentrations  were  reportedly  measured  with 
a  catalytic  type  hydrocarbon  sensor  and  recorded  at  2  Hz. 
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TABLE  IV. 3 

TEST  CONDITIONS  FOR  ESSO  TESTS  11  AND  17 


ESSO  11 

ESSO  17 

LNG  Volume  Spilled  (m3) 

10.2 

8.4 

Wind  Speed  @  5.6  m  (m/s) 

8.04 

4.1 

Atmospheric  Stability 

D 

D 

Ambient  Temperature  (K) 

300.2 

291.2 

Relative  Humidity  [%) 

78 

85 

Surface  Roughness  (m) 
(estimated) 

0.0001 

0.0001 

DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Esso  11. 


10  100  1000 

Distance,  m 

Figure  IV. 15.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Esso  17. 


IV. 4  U.S.  Department  of  Energy  Burro  and  Coyote  LNG  Releases 
on  Water 

The  U.S.  Department  of  Energy  sponsored  a  series  of  nine  LNG 
release  field  experiments  conducted  in  1980  at  China  Lake,  California 

3 

The  LNG  volumes  released  on  water  ranged  from  24  to  39  m  ,  at  rates 

3 

of  11.3  to  18.4  m  /min.,  with  subsequent  dispersion  over  land.  The 
extensive  data  from  the  experiments  were  reported  by  Koopman  et  al . 
(1982),  and  seyeral  analyses  of  the  test  data  and  comparisons  with 
predictions  using  dispersion  models  developed  by  Lawrence  Livermore 
National  Laboratories  (LLNL)  have  been  published  (Cederwell  et  al., 
1981;  Ermak  et  al . ,  1982;  Morgan  et  al.,  1983b;  Chan  and  Ermak,  1983) 

The  subsequent  Coyote  test  series,  designed  to  study  LNG  fire 
and  rapid  phase  transition  phenomena  also  provided  test  data  on 
atmospheric  dispersion  of  LNG  vapors.  The  test  data  for  the  Coyote 
series  are  reported  by  Goldwire  et  al .  (1983),  and  comparison  of 
the  experimental  vapor  dispersion  model  predictions  using  the  SLAB 
and  FEM3  models  developed  by  LLNL  have  been  published  (Morgan  et 
al.,  1983c). 

Table  IV. 4  shows  the  pertinent  test  descriptions  for  Burro  3, 

7,  8,  9,  and  Coyote  5  and  6,  which  have  been  simulated  with  the 
DEGADIS  model.  The  experiments  were  modeled  as  steady,  time-limited 
releases,  i.e.  Burro  9  was  represented  as  a  steady  release  of  130 
kg/s  for  79  seconds.  The  liquid  source  diameter  was  determined 
using  an  LNG-on-water  boiling  rate  of  0.085  kg/m  s,  suggested  by 
the  analysis  of  the  Shell  Map! in  Sands  LNG  release  experiments 
(Blackmore  et  al.,  1982).  In  all  cases,  LNG  is  assumed  characterized 
by  the  properties  of  methane.  The  reported  measured  concentrations 
are  based  on  a  10  s  averaging  time;  the  lowest  sensor  position  was 
1  m.  Ninety  gas  sensors  were  used  of  which  45  were  1ST  (Inter¬ 
national  Sensor  Technology) detectors ,  33  were  an  infrared 


104 


absorption  detector  developed  by  LLNL,  and  12  were  MSA  (Mine 
Safety  Appliance)  catalytic  detectors.  Figures  IV. 16  through 
IV. 21  show  the  OEGADIS  model -predicted  maximum  ground  level 
centerline  gas  concentration  compared  with  the  measured  (10  s 
averaged)  maximum  concentrations  as  a  function  of  downwind  distance 
from  the  release  center  for  these  tests. 


TABLE  IV. 4 


DOE  CHINA  LAKE  LNG  SPILL  TEST  SCENARIO  DATA  SUMMARY 


Test. 

'Ho. 

LN6 

Coap. 

t  c, 

*  c2 

*  C3 

Liquid  Liquid 
Spill  Spill 
Voluae  Bote 
tr>  or /min 

Wind 

Speed, 

a/s 

9  1  ■ 

3  m 

8  ■ 

T«ap. 

•c 

9  2a 

Ataos. 

Pressure 

ATM 

Nonin 

Obukhov 

Length 

a 

Pasqulll 

Stability 

Surface 

Teap. 

•c 

Relative 

Hualdity 

Upwind 

Downwind 

5.) 

Burro  3 

Unknown 

34.0 

12.2 

5.6 

33.8 

0.95 

-9.06 

C 

310 

-/5.2 

5.8 

92.9 

7.8 

Burro  7 

5.8 

39.4 

13.6 

8.8 

33.7 

0.94 

-114 

0 

310 

6. 7/7.4 

1.2 

9.6 

96.5 

1.6 

Burro  8 

2.8 

28.4 

16.0 

1.9 

33.1 

0.94 

♦16.5 

E 

310 

4. 7/4. 5 

0.6 

2.4 

93.6 

5.3 

Burro  9 

5.2 

24.2 

18.4 

5.9 

35.4 

0.94 

-140 

C 

310 

11.7/14.4 

0.9 

6.5 

74.9 

10.2 

Coyote  S 

20.5 

28.0 

17.1 

10.5 

27.9 

0.94 

-26.5 

C 

310 

22.1/  — 

4.6 

" 

11.9 

81 .8 

5.0 

Coyote  6 

14.6 

22.8 

16.6 

5.0 

24.4 

0.94 

♦73.6 

0 

310 

22.8/  -- 

3.6 

5.9 

100  1000 
Distance,  m 

DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Burro  3. 
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Figure  IV. 17.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Burro  7. 
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DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Burro  9. 


DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Coyote  5. 


Figure  IV. 21.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Coyote  6. 


IV. 5  Shell  Research  LNG  and  LPG  Releases  on  Water  at  Map! in  Sands 

Shell  Research  conducted  a  series  of  LNG  and  LPG  releases  on 
water  at  Maplin  Sands  in  1980-81  in  south  England.  Ten  continuous 
LNG  releases,  eleven  continuous  LPG  releases,  and  two 
instantaneous  releases  each  of  LNG  and  LPG  were  reported.  Some  of 
the  releases  were  ignited  to  obtain  cloud  and  pool  fire  heat  flux 
measurements,  whereas  some  experiments  were  conducted  primarily  to 
measure  gas  concentrations.  The  volumes  of  LNG  and  LPG  released 
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ranged  from  about  6  m  to  22  m  and  8  m  to  27  m  ,  respectively. 

Maximum  rates  of  release  in  the  continuous  LNG  and  LPG  spills  were 
3  1 

about  4.8  m  /min  and  4.3  m  /min,  respectively.  Data  from  these 
releases  have  been  reported  in  several  papers  by  Shell  personnel 
(Blackmore  et  al . ,  1981;  Puttock  et  al.,  1982a,  1982b,  1982c; 
Colenbrander  and  Puttock,  1983a),  and  extensive  data  digests 
summarizing  the  data  have  been  prepared  by  Shell  (Colenbrander  et 
al.,  1984).  Comparisons  of  the  experimental  data  with  predictions 
obtained  with  Shell's  HEGADAS  model,  which  was  the  starting  point 
for  the  DEGADIS  model  herein  proposed,  have  also  been  reported 
(Colenbrander  and  Puttock,  1983a). 

Table  IV. 5  shows  the  test  descriptions  for  six  continuous  LNG 
releases,  one  instantaneous  LNG  release,  and  six  continuous  LPG 
releases  conducted  at  Maplin  Sands  which  have  been  simulated  with 
the  OEGADIS  model.  The  continuous  releases  have  been  simulated 
as  steady  state  with  the  liquid-on-water-pool  size  determined 
assuming  boiling  rates  for  LNG  and  LPG  of  0.085  and  0.12  kg/m2s, 
respectively.  The  simulation  of  Maplin  22  assumes  a  steady, 
time-limited  release  of  510  kg/s  for  10  seconds.  All  simulations 
assumed  neutral  stability  with  a  reported’ surface  roughness  of 
3.38  x  10"4  m.  Figures  IV. 22  through  IV. 34  show  the  DEGADIS  model- 
predicted  maximum  ground  level  centerline  gas  concentration  compared 


TABLE  IV. 5 

SHELL  MAPLIN  SANDS  LN6/LPG  SPILL  TEST  SCENARIO  DATA  SUMMARY 


Test  No. 

Liquid 

Coinp. 

C1  C2 
c3  c4 

Liquid 

Spill 

Volume 

Liquid 

Spill 

Rate 

■Vain. 

Hind 
Speed 
•  10  a 

a/s 

Air 

Teap. 

8  2a 

#C 

Relative 

Huaidlty 

9  10  a 
t 

Surface 

Teap. 

#C 

Lowest 

Gas  Concentration 
Sensor 
a 

27 

93.2  5.4 
1.1  0.3 

12.6 

3.2 

5.5 

14.9 

53 

15.6 

0.9  -  1.0 

29 

98.5  1.4 
0.1  - 

21.9 

4.1 

7.4 

16.1 

52 

16.8 

0.9 

2 

-4 

M 

34 

95.9  2.6 
0.9  0.6 

10.2 

3.0 

8.6 

15.2 

72 

15.8 

0.9  >  1.0 

4* 

35 

97.8  1.7 
0.4  0.1 

18.3 

3.9 

9.8 

16.1 

63 

16.6 

0.9  -  1.0 

5 

39 

95.2  1.7 
0.6  0.2 

10.9 

4.7 

4.1 

16.7 

63*1  > 

— 

0.9  -  1.3 

56 

93.3  4.3 
1.5  0.4 

6.3 

2.5 

5.1 

10.6 

830) 

11.6 

0.6  -  0.7 

2 

—I 

22 

92.3  4.0 
0.9  0.2 

12 

WA 

S.l 

18.9 

62 

18.6 

0.9  -  1.3 

43 

—  1.22 
97.0  1.2 

17.2 

2.3 

5.5 

17.0 

~ 

18.9 

0.9 

46 

0.01  1.05 
97.3  1.6 

22.2 

2.8 

8.1 

18.7 

71 0) 

17.3 

0.9  -  1.0 

b 

J 

» 

| 

47 

0.01  1.78 
96.9  1.3 

17.6 

3.9 

5.6 

17.4 

78(1) 

17.1 

0.9  -  1.0 

1 

- 

49 

—  2.4 
96.5  0.66 

8.4 

2.0 

6.2 

13.3 

88 

13.0 

0.9 

5 

50 

-  1.18 
97.6  1.2 

17.2 

4.3 

7.9 

10.4 

79 

9.9 

0.9  -  1.3 

54 

-  1.1 
97.4  1.5 

11.7 

2.3 

3.8 

8.4 

«5»> 

9.4 

0.5  -  0.7 

(1)  Hualdlty  from  Foulness  Meteorological  Station  (5  ka  my  and  1  ka  Inland)  corrected  to 
Naplln  site  ta^erature. 
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with  measured  maximum  steady  state  concentrations  as  a  function  of 
downwind  distance  from  the  release  center.  All  reported  measured 
concentrations  are  based  on  a  3  s  moving  average.  The  aspirated 
concentration  sensor  used  a  balanced  Wheatstone  bridge  to  measure 
the  heat  loss  from  a  sensing  element  placed  in  the  sample  stream. 
Shaded  experimental  points  represent  maximum  concentrations  which 
occurred  before  or  following  the  time  interval  judged  to  best  delimit 
the  steady  release  condition. 
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Figure  IV. 22.  DEGADIS-predlcted  centerline  maximum  concentration  vs.  maximum  measured  concentratlon- 
Maplin  27. 
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.23.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
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Figure  IV. 24.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Maplin  34. 
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.25.  OEGADIS-predlcted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Maplin  35. 
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.26.  OEGAOIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
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Figure  IV. 28.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Maplin  22. 
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Figure  IV. 30.  DEGADIS-Jifedicted  centerline  maximum  concentration  vs.  maximum  measured  concentrSTtfon 
Hapltn46. 


OEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Maplln  47. 
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.32.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
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.33.  DEGADIS-predlcted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
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DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
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IV  J5  British -Health  and  Safety  Executive  Thorney  Island 
Heavy  Gas  Trials 

An  industry-government  consortium  organized  by  the  British 
Health  and  Safety  Executive  conducted  a  series  of  instantaneous 
releases  of  nominal  2000  m^  volumes  of  Freon-1 2/nitrogen  mixtures 
on  the  disused  Thorney  Island  airfield  on  the  southern  coast  of 
England  from  July  1982  -  July  1983  (McQuaid,  1983).  The  releases 
were  conducted  in  wind  speeds  (at  10  m)  ranging  from  1.7  m/s 
to  7.5  m/s  and  under  different  atmospheric  stability  conditions. 

In  all  of  the  experiments  the  gas  was  contained  in  a  twelve-sided, 
plastic-sheet-walled  bag  14  m  across  and  13  m  high.  The  bag  was 
supported  on  taut  rigging  from  a  central  support  column  with 
radial  guy  wires  to  ground  anchors  and  was  fitted  with  a  conical 
roof,  also  fabricated  from  plastic  sheeting,  which  was  withdrawn 
upwards  immediately  prior  to  release.  The  bag  was  filled  by 
vaporizing  liquid  Freon-12  and  liquid  nitrogen  and  heating  the  gas 
mixture  to  near  ambient  temperature.  The  releases  were  essentially 
isothermal.  The  sides  of  the  container  collapsed  in  concertina 
fashion  to  ground  level  in  less  than  two  seconds,  leaving  an 
upright  freestanding  cylinder  of  gas.  Detailed  descriptions  of  the 
experimental  design  have  been  published  as  have  complete  data 
digests  for  all  of  the  trials  (HSE,  1982-83).  Two  gas  sensor 
types  were  used;  the  standard  gas  sensor  was  an  oxygen  deficiency 
meter  based  on  a  semipermeable  membrane  and  electrochemical  cell. 

A  fast  response  sensor  (10  Hz)  basically  aspirated  a  sample  past 
the  cell  membrane.  Reported  measured  concentrations  are  based 
on  a  0.6  s  time  average.  The  lowest  sensors  were  0.4  m  above 
ground  while  the  grass  was  approximately  0.1  m  high. 


Table  IV. 6  shows  the  pertinent  test  descriptions  for  Thorney 
Trials  7,  8,  9,  11,  13,  and  15,  which  have  been  simulated  with  the 
DEGADIS  model.  Figures  IV. 35  through  IV. 40  show  thd  DEGADIS- 
predicted  average  source  concentration  for  shorter  distances  and 
the  maximum  centerline  gas  concentration  for  longer  distances 
compared  with  the  maximum  measured  concentration  as  a  function  of 
downwind  distance  from  the  release  center  for  these  trials. 


TABLE  IV. 6 

BHSE  THORNEY  ISLAND  PHASE  I  TRIALS— DATA  SUMMARY 


Test 

No. 

Volume 

Released 

(m3) 

Initial 

Relative 

Density 

o/Pa 

Wind 

Speed 

8  10  m 
m/s 

Air 

Temp. 

@  9  m 
°C 

Relative 

Humidity 

8 

10  m 

Pasquill 

Stability 

7 

2000 

1.75 

3.2 

17.1 

80.7 

E 

8 

2000 

1.63 

2.4 

17.1 

87.6 

D 

9 

2000 

1.60 

1.7 

18.6 

87.3 

E 

11 

2100 

1.96 

5.1 

12.3 

77.1 

D 

13 

1950 

2.00 

7.5 

13.2 

74.1 

D 

15 

2100 

1.41 

5.4 

10.3 

88.4 

C/D 

p 


w 

1 1 
3  ’« 


Distance,  m 

Figure  IV. 35.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Thorney  07 . 


Figure  IV. 36.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration 
Thomev  08. 
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Figure  IV. 37.  OEGAOIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentratlon- 
Thorney  09. 
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Figure  IV. 40.  DEGADIS-predicted  centerline  maximum  concentration  vs.  maximum  measured  concentration- 
Thorney  15. 
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IV. 7  Discussion  of  DEGADIS  Simulation  vs.  Data  Comparisons 


The  thirty-nine  field  experiments  simulated  include  releases 
exemplifying  a  wide  range  of  heavy  gas  dispersion  behavior.  The 
small  scale  propane  releases  reported  by  Welker  have  been  shown 
to  be  reasonably  well  described  by  a  passive  dispersion  model 
(Welker,  1982)  and  were  included  to  demonstrate  the  applicability 
of  DEGADIS  to  passive  dispersion-dominated  heavy  gas  releases. 

The  Thorney  Island  instantaneous  heavy  gas  releases  are  indicated, 
by  comparison  with  laboratory  calm  air  heavy  gas  releases  (Volume 
II,  this  report),  to  be  dominated  by  buoyancy-driven  flow  during 
the  period  when  the  cloud  concentration  decreases  to  about  5%.  The 
instantaneous  releases  onto  water  of  LNG,  as  well  as  the  quasi - 
continuous  releases  of  LNG  and  LPG  on  water,  reflect  important 
components  of  all  three  flow  regimes:  buoyancy-dominated  flow, 
stably  stratified  shear  flow,  and  passive  dispersion. 

Comparison  of  the  model  predictions  with  experimental 
observations  is  facilitated  by  classification  of  the  tests  with  a 
Release  Richardson  Number  which  represents  the  ratio  of  a 
buoyancy-driven  velocity  scale  characteristic  of  the  release  to  a 
turbulence  velocity  scale  characteristic  of  the  atmospheric  flow. 
Richardson  numbers  are  defined,  for  this  purpose,  for  continuous 
and  instantaneous  releases  as  follows. 


c  2 

Continuous  Releases:  Rig  =  gA'Q/uu*D 


I  1/3  2 

Instantaneous  Releases:  Rig  =  gA'  V/  /  u* 


where  the  characteristic  buoyancy-driven  velocities  are 


77T7T 


g  A'  H  =  /  g  A'  Q  /  u  D  and  /gA'  H  s  /gA1  V. 


for  continuous  and  instantaneous  releases  respectively,  and  the 
friction  velocity  is  assumed  characteristic  of  the  atmospheric 
flow. 


«\»\ 


V.'jYiViWW 


r 

Table  IV. 7  contains  a  summary  of  Rig  for  the  continuous  releases 
considered  as  well  as  an  indication  of  whether  a  secondary  gas 
"blanket"  was  formed  over  the  primary  gas  source  which  occurs  when 
the  atmospheric  takeup  rate  is  less  than  the  gas  evolution  rate 
(i.e.  when  the  atmospheric  shear  is  not  sufficient  to  assimilate 
the  gas  as  fast  as  it  is  being  produced  by  the  primary  source).  A 
secondary  gas  blanket  was  predicted  for  all  of  the  continuous 
releases  simulated  with  Richardson  numbers  above  32.  In  contrast 
to  the  high  dilution  rates  due  to  frontal  air  entrainment  in  the 
gas  blanket  for  releases  such  as  the  Thorney  Island  trials  and 
this  work  (Volume  II),  the  secondary  blankets  formed  for  these 
continuous  releases  are  not  predicted  to  be  significantly  diluted 

r 

by  air.  In  contrast,  the  low  Richardson  number  (Rig  =  0(1)) 
releases  are  predicted  to  have  sufficient  entrainment  of  gas  by 
the  wind  shear  to  not  only  preclude  the  formation  of  a  blanket,  but 
to  significantly  dilute  the  gas  in  the  shear  layer  over  the  source 
in  agreement  with  the  low  Richardson  number  range  indicated  above. 
Between  these  Richardson  number  extremes,  the  source  concentration 
is  predicted  to  decrease  with  decreasing  Richardson  number  without 
any  upwind  spreading. 

A  summary  of  the  instantaneous  release  Richardson  numbers  is 
given  in  Table  IV. 8  for  the  Thorney  Island  trials.  When  scaled  to  the 
Thorney  Island  releases,  the  laboratory  calm  air  releases  presented 
in  Volume  II  are  in  good  agreement  with  the  maximum  observed 
concentrations  for  the  Thorney  Island  trials  before  the  time  when 
the  average  concentration  has  dropped  to  about  5%.  This  suggests 
the  flow  and  dilutions  processes  during  this  period  reflect  a 
dominance  of  buoyancy-driven  flow  for  all  six  releases  (i.e.  for 
RiJ  >  0  (500)).  . 

In  order  to  quantify  the  uncertainty  associated  with  a 
predicted  distance  to  a  given  concentration  level,  a  quantitative 
comparison  between  the  predicted  and  "observed"  distance  to  the 
5%,  2-1/2%,  and  1%  concentration  levels  was  made.  The  "observed" 
values  were  determined  from  reported  maximum  concentrations  for 
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TABLE  IV. 7 

SUMMARY  OF  RICHARDSON  NUMBERS  FOR  CONTINUOUS  RELEASES 


Test 

o  o 

Blanket 

Formation 

Predicted 
Concentration 
at  Downwind 
Source  Edge 
(*> 

Welker  275-1 

0.1 

no 

1.8 

281-2 

0.6 

no 

296-1 

0.3 

no 

297-1 

0.2 

no 

297-3 

0.4 

no 

1.7 

298-1 

2.1 

no 

292-2 

1.3 

no 

2.4 

302-1 

0.1 

no 

334-1 

2.5 

no 

3.2 

334-2 

3.5 

no 

Burro  3 

29 

no 

83 

7 

8.9 

no 

55 

8 

930 

yes 

9 

32 

yes 

Coyote  5 

4.7 

no 

41 

6 

47 

yes 

Maplin  Sands  27 

18 

no 

77 

29 

9 

no 

57 

34 

5 

no 

46 

35 

4 

no 

42 

39 

55 

yes 

43 

34 

yes 

46 

12 

no 

43 

47 

42 

yes 

49 

21 

no 

63 

50 

15 

no 

46 

54 

100 

yes 

56 

20 

no 

80 

! 
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TABLE  IV. 8 

SUMMARY  OF  RICHARDSON  NUMBERS  FOR 
INSTANTANEOUS  RELEASES 


Test 


Thorney  Island  7 

8 

9 

11 

13 

850 

15 

690 

each  experiment  by  drawing  a  visual  best-fit  straight  line  (on  the 
Figures  in  this  Section)  through  the  reported  points  in  the  con¬ 
centration  range  of  interest.  All  of  the  measurements  used  were 
made  at  heights  at  or  below  1  m.  The  predicted  distance  to  a 
given  concentration  level  was  based  on  the  ground  level  centerline 
concentration  calculated  by  DEGADIS;  for  the  concentrations  and 
conditions  of  interest,  the  predicted  concentration  level  is 
essentially  constant  for  heights  below  1  m.  Table  IV. 9 
summarizes  the  "observed"  and  predicted  distances  to  the  5%, 

2-1/2%,  and  1%  concentration  levels  for  the  Burro,  Coyote,  Maplin 
Sands,  and  Thorney  Island  releases. 

Table  IV. 9  also  includes  ratios  of  the  "observed"  to  predicted 
distances  for  the  concentration  levels  of  interest.  As  well,  a 
90%  confidence  interval  (Conover,  1980)  of  these  ratios  is 
included  for  each  test  series  and  for  all  of  the  experiments 
together.  For  example,  the  Maplin  Sands  comparisons  indicate  the 
ratio  of  the  "observed"  to  the  predicted  distance  to  the  2-1/2% 
concentration  level  would  between  0.91  and  1.20  in  nine  out  of  ten 
realizations.  For  all  of  the  comparisons,  the  90%  confidence 
interval  of  the  ratio  of  the  "observed"  to  the  predicted  distance 
to  the  5%  level  would  be  0.73 and  0.96;  for  the  2-1/2%  level,  it 
would  be  0.82 and  1.03;  and  for  the  1%  level,  it  would  be  0.95  and 
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TABLE  IV. 9 

COMPARISON  BETWEEN  "OBSERVED"  AND  DEGADIS-PREDICTED  MAXIMUM  DISTANCE 
TO  GAS  CONCENTRATIONS  IN  THE  FLAMMABLE  CONCENTRATION  RANGE 


Distance  to  5%  Concentration 

Distance  to  2-1/2%  Concentration 

Distance  to  1%  Concentration 

(°es)J 

(PRE)f 

(08S/PRE)5 

(OBS)J  5 

(OBS/PRE)^ 

(OBS)f 

(PRE)f* 

(OBS/PRE), 

Test 

(■) 

(") 

(") 

(«) 

(■) 

(■) 

Burro  3 

200 

380 

0.52 

6 

320 

460 

0.69 

6 

540 

650 

0.83 

1 

7 

240 

380 

0.63 

2 

410 

550 

0.74 

5 

800 

840 

0.95 

2 

8 

420 

450 

0.93 

3 

740 

610 

1.21 

3 

1500 

900 

1.66 

7 

9 

240 

470 

0.51 

1 

480 

610 

0.78 

7 

1000 

820 

1.22 

0 

Coyote  5 

240 

320 

0.75 

0 

300 

420 

0.71 

4 

420 

590 

0.71 

2 

6 

210 

540 

0.38 

9 

300 

790 

0.38 

0 

470 

1080 

0.43 

5 

Burro  90% 
Confidence 

Interval 

0.46 

<  (0BS/PRE)s  _<  0.78 

0.55  <  (0BS/PRE)2 

5  <  0.98 

0.63  . 

c  (OBS/PRE) 

1  - 

31 

Maplln  22 

270 

570 

0.47 

4 

500 

770 

0.64 

9 

1200 

1040 

1.15 

4 

27 

200 

210 

0.95 

2 

320 

290 

1.10 

3 

600 

410 

1.46 

3 

29 

170 

190 

0.89 

5 

280 

290 

0.96 

6 

540 

450 

1.20 

0 

34 

180 

140 

1.28 

6 

290 

210 

1.38 

1 

600 

340 

1.76 

5 

35 

180 

140 

1.28 

6 

350 

220 

1.59 

1 

940 

360 

2.61 

1 

39 

140 

300 

0.46 

7 

260 

410 

0.63 

4 

680 

560 

1.21 

4 

43 

95 

130 

0.73 

1 

160 

190 

0.84 

2 

320 

350 

0.91 

4 

46 

120 

110 

1.09 

1 

200 

170 

1.17 

6 

400 

300 

1.33 

3 

47 

140 

180 

0.77 

8 

240 

260 

0.92 

3 

500 

470 

1.06 

4 

49 

140 

110 

1.27 

3 

230 

170 

1.35 

3 

440 

300 

1.46 

7 

50 

120 

130 

0.92 

3 

200 

220 

0.90 

9  ■ 

480 

380 

1.26 

3 

54 

180 

160 

1.12 

5 

280 

230 

1.21 

7 

490 

410 

1.19 

5 

56 

150 

160 

0.93 

8 

230 

250 

0.92 

0 

400 

390 

1.02 

6 

Maplln  90% 
Confidence 

Interval 

0.80 

<  (obs/pre)5  <  1 

.10 

0.91  < 

(0BS/PRE)2  5  <  1. 

20 

1.15 

<  (OBS/PRE),  <  1 

.47 

Thorney  Island  7 

110 

94 

1.17 

0 

180 

160 

1.12 

5 

380 

440 

0.63 

6 

8 

120 

no 

1.09 

1 

190 

210 

0.90 

5 

390 

590 

0.66 

1 

9 

130 

130* 

1.00 

0 

190 

200+ 

0.95 

0 

320 

250 

1.28 

0 

11 

95 

170 

0.55 

9 

160 

260 

0.61 

5 

300 

470 

0.63 

8 

13 

120 

170 

0.70 

6 

190 

270 

0.70 

4 

340 

480  . 

0.70 

8 

15 

190 

180 

1.05 

6 

270 

300 

0.90 

0 

440 

570 

0.77 

2 

Thorney  Island  90% 
Confidence  Interval 

0.71 

<  (08S/PRE)5  <  1 

.11 

0.70  <  (08S/PRE)2 

.5  —  1 ' 

02 

0.64 

<  (OBS/PRE),  <  0.99 

Suamary  90% 
Confidence 

Interval 

0.73 

<  (0B$/PRE)5  <  0.96 

0.82  < 

(obs/pre)2>5  <  1. 

03 

0.95 

<  (OBS/PRE),  <  1 

.24 

Suamary  99% 
Confidence 

Interval 

0.64 

<  (obs/pre)5  <  1 

a 

BE 

(0BS/PRE)2  5  <  1. 

09 

0.88 

<  (OBS/PRE),  <  1 

.33 

‘"Observed"  distances  to  the  given  concentration  level  are  based  on  a  visual  best-fit  straight  line  through  the 
maxim*  reported  concentrations  for  each  of  the  trials. 

“OESADIS -predicted  distances  to  the  given  concentration  level  are  based  on  the  ground  level  centerline 
concentrations. 

*0E6A0IS-pred1cted  distance  to  the  given  concentration  level  based  on  twice  the  average  concentration  of  the 
"source"  cloud. 


1.24.  As  an  example,  if  for  a  given  release  scenario,  the  predicted 
distance  to  the  2-1/2%  concentration  level  was  120  m,  the  distance 
to  the  2-1/2%  concentration  level  for  nine  out  of  ten  realizations 
would  be  expected  to  range  between  98  m  and  124  m. 

In  all  of  these  simulations,  DEGADIS  used  correlations  for 
gravity  spreading  and  frontal  entrainment  during  the  buoyancy- 
dominated  flow  regime  based  on  the  laboratory  data  presented  in 
Volume  II.  For  stably  stratified  shear  flow,  the  mixing  data  of 
McQuaid  (1976),  Kantha  et  al.  (1977),  and  Lofquist  (1960)  were  used 
as  discussed  in  Section  I.  Heat  transfer  from  the  substrate  to  the 
cloud  was  based  on  correlations  for  heat  transfer  from  a  flat 
plate  as  developed  in  Section  I II . 1 .  Although  field  data  have  been 
used  to  describe  the  crosswind  passive  dispersion  parameter  ay 
(Pasquill,  1983)  and  the  along-wind  passive  dispersion  parameter  ax 
(Beals,  1971),  the  model  applicability  is  better  justified  if  no 
adjustment  to  laboratory  correlations  is  necessary  for  simulations 
to  be  consistent  with  observed  field  data.  Accordingly,  no  such 
adjustments  have  been  made. 

IV. 8  Effect  of  Parameter  Variation  on  DEGADIS  Simulations 

The  steady  state  base  case  HTAG  dispersion  scenario  described 
in  Table  IV. 10  was  simulated  to  provide  an  illustration  of  the  effect 
of  variation  of  the  important  scenario  definition  parameters: 
type  of  gas,  wind  speed,  atmospheric  stability,  surface  roughness, 
ground- to-cloud  heat  transfer,  and  atmospheric  humidity. 

TABLE  IV. 10 

STEADY  STATE  BASE  CASE  DISPERSION  SCENARIO 

Type  of  Gas:  LNG  vapor  3 

Source  Rate:  127.5  kg/s  (0.3  nr/s  liquid) 

Atmospheric  Stability:  D  (Pasquill) 

Surface  Roughness:  10"^  m 
Ground-to-Cloud  Heat  Transfer:  DEGADIS 
correlation 

Air  Temperature/Surface  Temperature:  25°C 
Atmospheric  Humidity:  50% 


Table  IV. 11  shows  the  effect  of  variation  of  the  prescribed 
parameters  from  the  steady  state  base  case  on  the  predicted 
maximum  distance  to  the  5%  gas  concentration  level.  When  changing 
atmospheric  stability  class  and  surface  roughness,  several 
default  parameter  values  are  changed  including  the  Monin-Obukhov 
length  (which  affects  ct  and  u*)  as  well  as  the  along-wind  and 
crosswind  passive  dispersion  parameters.  Noting  that  the  friction 
velocity  is  a  measure  of  the  ambient  turbulent  kinetic  energy  and 
therefore  a  measure  of  the  ambient  energy  available  to  disperse 
the  contaminant,  the  downwind  LFL  distance  (5%  concentration 
level)  is  inversely  proportional  to  the  friction  velocity  for  the 
same  wind  speed  for  the  0  stability  cases  (i.e.  when  a  is  the 
same  value).  Therefore,  in  the  absence  of  obstacles,  an  increase 
in  surface  roughness  would  result  in  a  concomitant  decrease  in 
LFL  distance.  In  comparing  the  methods  of  estimating  the  heat 
transfer,  wind  speed  and  ambient  humidity  both  affect  the  result. 

At  high  wind  speeds  (10  m/s),  the  LFL  distance  is  relatively 
insensitive  to  how  the  heat  transfer  is  estimated  (or  if  heat 
transfer  is  included  at  all).  For  high  atmospheric  absolute 
humidities,  the  LFL  distance  is  insensitive  to  wind  speed  under 
these  conditions;  this  behavior  is  not  observed  for  propane. 

The  effect  of  initial  volume  and  wind  speed  on  the  maximum 
downwind  distance  to  the  5%  concentration  level  for  instantaneous 
LNG  releases  of  10  m3,  100  m3,  1,000  m3,  and  10,000  m3  is  shown  in 
Table  IV. 12  and  Figure  IV. 41.  For  the  release  conditions  simulated 
there  is  comparatively  little  difference  between  the  maximum  LFL 
distance  for  the  1.25  m/s  and  2.5  m/s  wind  speeds  for  the  initial 
volumes  simulated.  For  the  largest  initial  volume  simulated,  there 
is  also  little  difference  between  the. predicted  distance  to  the  5% 
concentration  level  for  the  2.5  m/s  and  5  m/s  wind  speeds. 


TABLE  IV. 11 

EFFECT  OF  SOME  INPUT  PARAMETERS  ON  THE  MAXIMUM  DOWNWIND  DISTANCE 
TO  THE  5%  GAS  CONCENTRATION  LEVEL 
FOR  A  STEADY  LNG  RELEASE 


Mind  Speed: 
2.5  m/s 
Distance  (m)‘ 

Mind  Speed: 

5  m/s 

Distance  (m)* 

Mind  Speed: 
10  m/s 

Olstance  (m)‘ 

Base  Case 

860 

600 

400 

Parameters 

Isothermal  Gas  , 

(p,  -1.79  kg/m  ) 

720 

430 

220 

Atmospheric  Stability 

Pasqulll  B 

500 

370 

230 

Pasqulll  F 

1380 

— 

— 

Surface  Roughness 

10-J  m 

700 

500 

320 

IQ'5  m 

940 

700 

490 

Substrate  Heat  Transfer 

None 

1370 

810 

470 

Using  LLNl  Correlation  400 

430 

400 

Atmospheric  Humidity 

80S 

350 

380 

300 

20X 

670 

620 

410 

‘Release  conditions  given  In  Table  IV. 10. 


TABLE  IV. 12 

EFFECT  OF  WIND  SPEED  AND  VOLUME  ON  THE 
PREDICTED  MAXIMUM  DOWNWIND  DISTANCE  TO  THE 
5%  GAS  CONCENTRATION  LEVEL  FOR  AN 
INSTANTANEOUS  LNG  RELEASE  ON  WATER 


Liquid 

Volume 

(m3) 

Mind  Speed: 
1.25  m/s 
Olstance  (m)* 

Mind  Speed: 

2.5  m/s 
Olstance  (m)* 

Mind  Speed: 

5  m/s 

Distance  (m)+ 

Mind  Speed: 
10  m/s 

Olstance  (m)‘ 

10 

640 

550 

480 

320 

100 

1330 

1350 

1120 

800 

1,000 

3100 

2830 

2220 

10,000 

7000 

6800 

5300 

‘Release  Conditions:  Pasqulll  0  stability,  air  temperature  of  298°K, 
68X  relative  humidity,  sea  teaeerature  of 
288°K,  LNG  belling  rate  of  0.065  kg/m*  s. 


*  LNG  Volume  (m°) 

Figure  IV. 41.  Effect  of  wind  speed  and  volume  on  the  predicted  maximum  downwind  distance 
to  the  5%  gas  concentration  level  for  an  Instantaneous  LNG  release  on 
water. 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


Based  on  review  of  the  phenomenology  of  heavier-than-air  gas 
(HTAG)  dispersion,  the  following  phases  are  used  to  describe  HTAG 
dispersion  in  the  atmosphere: 

—  Negative  buoyancy-dominated  dispersion 
—  Stably  stratified  shear  flow 
—  Passive  dispersion  due  to  atmospheric  turbulence 
Modeling  concepts  based  on  investigations  of  laboratory  gravity 
currents  have  been  used  to  describe  the  important  HTAG  flow  and 
dilution  processes  that  characterize  buoyancy -dominated  dispersion. 
Investigation  of  the  gravity  spreading  and  dilution  of  HTAG  volumes 
released  suddenly  in  calm  air  demonstrates  the  applicability  of 
proposed  scaling  methods;  also  for  these  releases,  no  significant 
change  in  measured  peak  concentration  as  a  function  of  distance 
was  observed  over  the  tested  range  of  initial  height-to-diameter 
ratio  (0.4  _<  (H/D) ^  <  1.57)  and  initial  specific  gravity 
(2.16  <  (p/pa)^  <  4.19).  Scaling  of  laboratory  experimental  data 
from  releases  of  small  volumes  (35  liters  to  530  liters)  to 
similar  releases  of  up  to  2000  m  in  the  Thomey  Island  Trials 
has  been  demonstrated. 

Laboratory  data  from  stratified  shear  flow  mixing  experiments 
have  been  used  to  model  the  vertical  dispersion  of  HTAG  in  the 
atmospheric  constant  stress  layer  consistently  with  the  limiting 
passive  dispersion  behavior  of  demonstrated  air  pollution  models. 

An  interactive  computer  model  for  HTAG  gas  dispersion 
(DEGADIS),  which  can  simulate  a  wide  variety  of  HTAG  release 
scenarios,  including  instantaneous  releases,  steady  releases,  and 
prescribed  time-varying  releases,  has  been  developed.  The  model 
accounts  for  the  three  regimes  of  dispersion  described  above  and 
for  effects  due  to  energy  exchange  between  the  dispersing  cloud 


and  the  underlying  surface.  DEGAOIS  estimates  an  ensemble-average 
concentration  downwind  of  a  ground  level  area  for  flat,  homogeneous 
terrain  (I.e.  a'maxlmum  surface  roughness  of  order  0.01  m). 

The  DEGAOIS  model  has  been  used  to  simulate  a  wide  range  of 
field  experimental  HTAG  releases.  Including  small  to  Intermediate 
LPG  (0.1  to  1  kg/s)  and  LNG  releases  (1  -  100  kg/s)  on  land,  large 
scale  releases  (10  -  150  kg/s)  of  LPG  and  LNG  on  water,  and 
instantaneous  releases  of  approximately  5000  kg  Freon/air  mixtures 
on  land.  The  DEGADIS  model  predicted  downwind  gas  concentration 
decay  has  been  shown  consistent  with  data  from  the  field  tests 
considered  representing  the  maximum  range  of  field  experimental 
data  currently  available.  The  DEGADIS  model -predicted  source 
characteristics  have  been  shown  consistent  with  the  continuous 
water  tunnel  experiments  of  Britter  (1980a)  based  on  a  release 
Richardson  number  reflecting  the  relative  importance  of  the 
negative  buoyancy-driven  flow  to  the  ambient  shear  flow. 

Applicability  of  DEGADIS  has  been  primarily  directed  to  the 
prediction  of  concentrations  In  the  lower  flammability  limit 
range  (-1  to  5*).  Even  though  the  relation  between  the  peak  and 
ensemble-averaged  gas  concentration  Is  uncertain,  there  Is  some 
basis  for  using  2.0  as  an  estimate  of  the  peak-to-ensemble- 
averaged  concentration  ratio  for  determining  a  flammable  gas  con¬ 
centration  zone;  the  maximum  distance  to  a  flammable  gas  concentra¬ 
tion  would  then  be  the  distance  predicted  to  one-half  the  lower 
f 1  ammabi 1 i ty  limit  (LFL/2).  Based  on  comparison  with  field  data, 
the  ratio  of  observed  to  calculated  distance  for  the  2-1/2% 
concentration  would  be  expected  to  range  from  0.78  to  1.09  for  99 
out  of  100  realizations  (I.e.  a  99%  confidence  interval);  the  99% 
confidence  Interval  for  this  ratio  was  between  0.64  and  1.03  for 
the  5%  level  and  between  0.88  and  1.33  for  the  1%  level.  If  for 
a  given  LNG  release  scenario,  the  calculated  distance  to  the  2.5% 
ensemble-average  concentration  level  was  120  m,  the  downwind 
extent  of  a  flammable  gas  concentration  would  be  expected  to  range 
between  94  m  and  131  m  for  99  out  of  100  realizations  if  the 
peak-to-ensemble-averaged  ratio  of  2  Is  assumed. 
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The  consistency  of  the  model  predictions,  which  reflect  small 
scale  laboratory  fluid  flow  and  mixing  data,  with  the  results  of  a 
wide  range  of  field  experimental  releases,  suggests  the 
applicability  of  the  model  to  prediction  of  dispersion  from  much 
larger  releases.  The  model  is  recommended  on  this  basis  for  use 
in  predicting  dispersion  from  very  large  releases  of  flammable 
gas,  extending  to  the  largest  single  container  volumes  (-25000  m3) 
presently  in  marine  use. 

The  model  can  readily  be  modified  to  incorporate  better 
turbulent  mixing  and  heat  transfer  sub-models.  It  is  likely  that 
improvements  in  these  two  areas,  which  may  result  from  additional 
experimental  data  as  well  as  continuing  analysis  of  the  large  data 
base  which  has  become  available  during  this  work,  can  be  expected 
to  provide  the  most  important  information  for  improving  the 
confidence  level  in  heavy  gas  dispersion  prediction. 

The  model  has  been  tested  primarily  for  the  prediction  of 
dispersion  of  heavy  gas  clouds  to  concentrations  characteristic  of 
the  lower  flammability  limits  of  hydrocarbon  gases,  and  has  not 
been  demonstrated  for  application  to  the  prediction  of  dispersion 
to  the  ppm  concentration  range  which  characterizes  the  lower  limits 
of  toxic  gas  hazard.  However,  the  model  has  been  developed  to  be 
consistent  with  the  existing  data  base  on  passive  turbulent 
dispersion  in  the  far  field,  and  it  is  recommended  that  the  model 
be  evaluated  against  experimental  data  which  have  **  available 
for  dispersion  of  ammonia  and  nitrogen  tetroxide. 
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APPENDIX  A 

CLOUD  ENTHALPY  AND  DENSITY  CALCULATIONS 

The  estimation  of  cloud  density  is  an  important  part  of  the 
DEGADIS  model,  especially  for  releases  of  a  cryogenic  gas  with 
molecular  weight  lower  than  that  of  air.  With  regard  to  the  density/ 
enthalpy  calculations,  DEGADIS  can  be  run  in  one  of  three  modes. 

If  no  heat  transfer  is  to  be  included  in  the  calculations,  the 
density  can  be  specified  as  a  function  of  concentration  in  two 
different  ways.  When  either  the  release  is  isothermal  (Thorney 
Island  Trials)  or  the  user  needs  to  take  into  account  some  chemical 
reaction  (nitrogen  tetroxide)  or  phase  change  (ammonia)  of  the 
contaminant,  the  density  as  a  function  of  concentration  is  input  as 
a  list  of  data  triplets  consisting  of  mole  fraction,  concentration 
(kg/m  ),  and  mixture  density  (kg/m  ).  Linear  interpolation  is  used 
between  points. 

For  cryogenic  gas  simulations  without  heat  transfer,  the 
enthalpy  and  density  relationships  are  determined  by  adiabatic 
mixing.  For  ideal  solutions  of  ideal  gases  with  the  same  constant 
molal  heat  capacity,  adiabatic  mixing  gives 


Y  =  constant 


(A-l) 


where  p  is  the  mixture  density  and  cc  is  the  concentration  of 
contaminant.  Since  any  possible  condensation  of  ambient  humidity  is 
not  included  in  Equation  (A»l),  a  more  general  procedure  is  used. 

For  ideal  solutions  of  mixtures  of  contaminant/air/water,  the  mixture 
enthalpy  h(T)  can  be  expressed  by 

h(T)  =  wc  Cp  (T)  (T  -  Tr)  +  wA  (1  +  Hum  (T) )  C  (T  -  TR) 
c  Ka 

+  ("w  -  wA  Hum  (T))  [C  (T  -  TR)  -  X1]  (A-2) 

pw 


where  wc>  w^,  and  ww  are  the  mass  fractions  of  contaminant,  air,  and 
water,  respectively.  The  heat  capacities  of  air  and  water  are 
assumed  constant  at  1006  and  1865  J/kg°K,  respectively,  while  the 
heat  capacity  of  contaminant  is  described  (after  Hougen  et  al.,  1967) 
by  a  mean  molal  heat  capacity  given  by 


c  (T)  *  MWC  C  (T)  =  MW 
Kc,m 


3.33  x  10*  +  q.| 


r  P]  pi 

T  -  T 

1  'o 


T  -  T- 


(A- 3 ) 

where  Tg  is  the  contaminant  storage  temperature  (K)  and  values  of  q-j 
and  Pi  are  given  in  Table  A-l  for  the  range  90  K  _<  T  360  K.  Note 
that  a  constant  mean  molal  heat  capacity  can  be  used  by  setting 
Pi  =  1  and  qi  to  the  necessary  value.  The  latent  heat  given  in 
Equation  (A-2)  is  modeled  as  a  ramp  function  around  the  freezing 
point  of  water  by 


X' 


’  Avap 
Xvap+  Xfus 
3  Xvap  +  xfus 


T  -  263.151 

- TO - 

»  / 


for  T  273.15  K 

for  263.15  _<  T  <  273.15 

for  T  <  263.15  (A-4) 


In  order  to  determine  the  absolute  humidity  to  be  used  in  Equation 
(A-2),  the  saturation  conditions  are  first  determined.  The  saturation 
vapor  pressure  of  water  p*  is  estimated  with 

W 


p*  =  6.0298  x  10“3  exp 

W 


5407 


'  1 
273.15 


11 

T 


(A-5) 


where  p,*  is  in  atmospheres.  The  saturation  conditions  are  then  given 

W 


1 


A- 3 


where  p  is  the  total  pressure.  The  absolute  humidity  in  Equation 
(A-2)  is  the  minimum  of  the  saturated  humidity  and  the  ambient 
humidity,  i.e.  Hum  3  mi ri  (H$at,  Hfl).  Finally,  the  reference 
temperature  is  chosen  as  the  ambient  air  temperature. 

The  sum  of  the  pure  component  enthalpies  in  a  particular 
mixture  then  determines  the  enthalpy  of  the  mixture  under  the 
assumption  of  adiabatic  mixing.  Due  to  the  choice  of  reference 
temperature,  the  enthalpy  associated  with  any  entrained  air  is  zero. 

A  trial  and  error  procedure  calculates  the  mixing  temperature  from 
its  enthalpy. 

For  cryogenic  gas  simulations  with  heat  transfer,  the  amount  of 
heat  added  to  the  mixture  (D^)  is  calculated  by  the  model.  This 
determines  the  enthalpy  of  a  given  mixture  by 

h(T)  =  hAM  +  Dh  (A-7) 

where  hAM  is  the  adiabatic  mixing  enthalpy  of  a  given  mixture.  A 
similar  trial  and  error  procedure  then  determines  the  mixture 
temperature  T  from  the  new  value  of  enthalpy. 

With  the  mixture  temperature  known,  the  density  of  the  mixture 
is  calculated  assuming  Amagat's  law  of  additive  volumes  for  the  vapor 
phase  given  by 

p  =  [V  +  7  +  y  ]_1  (A-8) 

p  L  a  w, vapor  ycJ  v  ' 

where  V ,  +  V  is  the  volume  of  air  and  water  in  the  vapor  phase 

a  w, vapor 

given  by 


va  +  ^w, vapor  =  wa(°-00283  4  0.00456  Hum)  T/p 
The  volume  of  contaminant  in  the  vapor  phase  is 


V 


c 


w_ 


T 

PQ^O 


(A-9) 


(A- 10) 
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where  pg  is  the  density  of  the  saturated  vapor  at  the  storage 
temperature  Tg  and  ambient  pressure  p. 


TABLE  A-l 

VALUES  FOR  MEAN  MOLAL  HEAT  CAPACITY  CONSTANTS 


Species 

Pi 

Pi 

methane 

5.00 

5.6  x  10~8 

ethane 

2.79 

0.266 

propane 

2.25 

15.4 

